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ABSTRACT 


The problem of steady incompressible flow for lifting 

surfaces is considered. This problem requires the solution 

of an integral equation relating the values of the potential 

discontinuity on the lifting surface and its wake to the 

values of the normal derivative of the potential which are 

known from the boundary conditions. The lifting surface 

and the v?ake are divided into small quadrilateral (hyper- 

boloidal) surface elements, I., which are described in 

1 

terms of the Cartesian components of the four corner points. 
The values of the potential discontinuity and the normal 
derivative of the potential are assumed to be constant 
within each lifting surface element and equal to their 
values at the centroids of the lifting surface elements. 

This yields a set of linear algebraic equations. 

An iteration procedure is used to obtain the v?ake 
geometry: the velocities at the corner points of the wake 

elements are calculated and the (originally straight) wake 
streamlines are aligned -to be parallel to the velocity vec- 
tor . The procedure is repeated until convergence is attained 
Niomerical results are in reasonable agreement with exist 


rng ones. 
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SECTION I 


INTRODUCTION 


1.1 Definition of the Problem 


This work deals v/ith a nonlinear finite-element analy- 
sis of zero thickness wings (lifting surfaces) in steady, 
incompressible, inviscid, irrotational flow, including the 
effect of the rolled -up wake. The problem is formulated in 
terms of the velocity potential. This problem was consi- 
dered in Ref. 1 where a zeroth order finite-element (i.e. 
the potential (Ji prescribed at the centroids of the surface 
elements) analysis was used, with a straight- vortex- line 
wake. The present work is an extension of Ref. 1 and in- 
cludes the analysis of the wake roll-up as well as the non- 
linearities in the evaluation of the pressure (Bernoulli's 
Theorem) . Throughout this work, the potential is assumed 
to have a constant value over the surface element, equal for 
example with its (unknown) value at the centroid of the 
element (zeroth order formulation) . The first item considered 
here is the wake roll-up. The rolled-up geometry for the 
wake is obtained by a process of iteration. The convergence 
of the iteration scheme is investigated. A second item 
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included here is the effect of the rolled-up wake on the 
pressure distribution over the lifting-surface, using the 
nonlinearized Bernoulli equation. The results are compared 
with the linearized ones . 

1. 2 Lifting-Surface Theory 

The theoretical investigation of pressure and lift 
distributions over lifting surfaces of various shapes is 
embodied in many works.’ An excellent review of the litera- 
ture in the field is given in Refs. 2 and 3, together with 
results for lifting surfaces in steady and oscillatory, 
subsonic and supersonic flows. It may be worth noting that 
the integral equation used here is analogous to the one used 
by Jones {Ref. 4) for .unsteady incompressible flow. The 
classical approach for the n\unerical solutions of lifting- 
surface theories is by expressing the unknown in terms of a 
series with N unknown coefficients and by imposing that the 
equation be satisfied at N control points. Recently, however, 
a new approach (often referred to as the finite-element me- 
thod) has been introduced, especially in .connection with 
complex-configurations aerodynamics. A finite-element ana- 
lysis of lifting surfaces is considered for instance in Ref. 

5, which presents results for the loading of a rectangular 
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planar aero-thickness wing using a dovrnwash-velocity poten- 
tial formulation. Ref. 1 presents a general finite-element 
solution of a velocity potential formulation for lifting 
surfaces of arbitrary shapes in steady subsonic flows.. This 
work differs from the one of Ref. 5 in that it uses hyper- 
boloidal (i.e. warped) quadrilateral elements and is there- 
fore applicable to any arbitrary nonplanar shape. Expressions 
for the velocity at any point are also obtained in Ref.' 1. 
These are suitable for investigating the dynamics of the 
wake. 

1.3 Wake Roll-Up 


The interest in the phenomenon of wake roll-up has been 
spurred by the introduction, a few years ago, of the wide- 
body aircraft. Many papers have since been written about 
wing-tip vortices: about their formation, their effect on 
a trailing aircraft, their detection and their disappearance. 
Excellent descriptions of the phenomenon can be found in Refs. 
6, and 7. A short illustration of it is also presented here. 

In a few words, behind every aircraft in flight, a pair of 
counterrotating wing-*tip vortices is formed. See Fig. 1. 

The diameter of the vortex core has been found by measure- 
ments to be approximately 3% of the wingspan. The strength 
of the vortex seems to increase as the weight of the air- 
craft increases. If a four-engine j'et airplane flies 
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sufficiently high for the contrails to appear, it is ob- 
served that the exhausts from the two engines on each wing 
are gradually pushed towards the v;ing tips, thus making 
the wing-tip vortices visible. These vortices are quite 
stable; vortex life spans of more than 15 minutes have 
been observed, which, compared with the speed of a modern 
aircraft, means that the wing-tip vortices might persist 
for 150 miles behind the generating aircraft. The circum- 
ferential velocity of the vortex is large, of the order of 
30% of the generating aircraft speed. If a small aircraft 
passes through the wake of a large one, structural damage 
may occur on the small plane; if the flight path is not suf- 
ficiently high, the disturbances induced by the wake of the 
large aircraft on the velocity field of the small one may 
lead to loss of lift for the small plane and possibly to 
its crash. Ref. 7 contains more descriptive and pictorial 
information about these undesirable occurrences. 

Numerous wind tunnel and real life measurements of the 
wake vortices have been performed. See, for example, Ref s . 8 
and 9. 

Theories dealing with the matter are mainly two-dimen- 
sional and generally they do not account for the viscosity 
effects {Ref. 10). Ref. 11 presents a three-dimensional 
potential method for the estimation of the wake roll-up 
geometry for wings with control surfaces. In addition, an 
"artificial" viscosity coefficient is introduced in the 
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eguations describing the velocity field of the vortex sheet 
to "smoothen" out the singularities inherent in the method. 

Reference 12 presents another three-dimensional potential 
model to obtain a rolled-up wake geometry, as v;ell as the 
wing-jet interaction. 

Reference 13 integrates in time a set of ordinary differ- 
ential equations describing the position of the wake vortices. 
The finite vortex sheet of which the. wake consists is approxi- 
mated by a finite number of vortices. An unsatisfactory wake 
pattern was obtained and the paper contends that the mathema- 
tical model used fails at the wing tip. 

Reference 14 presents a method for the prediction of the 
aerodynamic loads on thin lifting surfaces. Nonlinearities 
(wake deformation) are considered. The method of Reference 
14 is conceptually the closest to the one presented in this 
work. 

1. 4 Formulation of the Problem 

This subsection presents the, basic flow equations which 
will be used throughout the paper. The fluid considered here 
is incompressible, inviscid and irrotational. For an in- 
compressible fluid, the continuity equation is 


? • if = 0 


( 1 . 1 ) 




- 6 - 


where V is the velocity vector. Because of the fact that 
the fluid is irrotational , or 


V X V = 0 (1.2) 

a velocity potential f exists, such that 

V = V$ (1.3) 


It is convenient to introduce the perturbation velocity 
potential ,<}> , and def ine V as 

V = (1 +V(f)) (1.5) 


where x is the unit vector along the x-direction. 

Combining now Eqs . • (1.1) and (1.5) the Laplace equation for ^ 
is obtained: 


V^(|) = 0 (1.6) 

The boundary condition to be satistfied is that the flow 
is tangent to the surface, or 

V • n = 0 (1.7) 

From Eqs. (1.5) and (1,7), the boundary condition for the 
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perturbation potential results : 

( 1 . 8 ) 

As it is well known, on the surfaces of the wing and of 
the wake the solution is discontinuous (see for instance 
Ref. 1 and 5). Also, there exists a pressure discontinuity 
on the surface of the wing, while the surface of the wake is 
determined by the fact that no pressure, discontinuity exists 
on the wake. Therefore, in order to complete the problem, 
the condition for the geometry of the wake as well as the 
expressions for the pressure discontinuity on the wing are 
obtained here. This can be easily accomplished, starting 
from the Bernoulli theorem (for steady, incompressible, invis- 
cid flows) 



P “ Pco 



(V 



(1.9) 


If there exists a surface of discontinuity, then, indi- 
cating for simplicity with "upper" and "lower" the two sides 
of the surfaces, one obtains, from Eq. (1.9) 


■u 



V - V • V ) = 
u a 


0 


( 1 . 10 ) 
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or 


Pu *- 


+ f (V. 


u 




(V - 
u 


Vj) . 0 


( 1 . 11 ) 


• J/iddcating with the velocity of the point on the surface 
of discontinuity 



V 

u 


+ V, 


( 1 . 12 ) 


(average between the upper and the lower surface) and with 




(1.13) 


the velocity discontinuity, Eq. (1.11) may be rewritten as 

Ap = - P^= -pV^ • AV (1.14) 


This is the desired expression for the pressure discontinuity. 
Using Eq. (1.5), Eq. (1.14) may be rewritten as 


Ap = -pu;; 


(i 




V(A<|)) 


(1.15) 


or 
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-2(1 + ^(f) ) . V(A4>) 

a 


(1.16) 


which gives the exact (nonlinear) expression for the pressure 
distribution on the wing. 

Equation (1.14) may also be used to obtain the condi- 
tion for the geometry of the wake. For, the condition that 
no pressure discontinuity exists on the wake yields 

V • Ay - 0 (1.17) 

a - 

It may be noted that if Eq. (1.17) is satisfied, then 
the no-pressure-discontinuity condition is automatically satis- 
fied. Equation (1.17) may be interpreted as saying that the 
velocity discontinuity on the wake is normal to the velocity 
of the wake. Also, Eq. (1.17) may be rewritten as 

(V • V) A(f) = 0 (1.18) 

d. 


i . e . that 


A(j) = constant along a streamline (1.19) 

Therefore, the geometry of the wake may be obtained from the 
streamlines emanating from the trailing edge which have the 
property of being tangent to V. Equation (1.17) (and hence 
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the condition of no-pressure-discontinuity) is then satisfied 
by imposing that A<J> be constant along a streamline (Eg. (1.19)) . 
It may be worth mentioning that Eg. (1.17) is- eguivalent to 
saying that the vortex lines coincide with the streamlines 
since a surface of velocity discontinuity (with continuous 
normal component) is eguivalent to a layer of vortices with 
vortex lines parallel to the lines of constant (which, 
in turn, are normal to the directions of AV)*. 

It is worth noting that the above formulation is exact, 
in the sense that no small-perturbation hypothesis has been 
used. In order to assess the relevance of using the exact 
formulation, the results obtained with such a formulation 
will be compared with the ones obtained from a small-pertur- 
bation formulation. If the small-perturbation hypothesis 

= 0(e) << 1 (1.20-) 

is invoked. Eg. (1.5) yields 

V = U 1 + 0 (e.) (1.21) 

<x> 

and therefore Eg. (1.16) may be rewritten as' 

*See for instance Ref. 15 
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ACp =:-2i • V (A(()) + 0{e) 


( 1 . 22 ) 


while the wake may be assumed to be composed of straight 
vortex lines emanating from the trailing edge. A more con- 
venient expression for ACp is 


+ 0(E) 


(1.23) 


where s is the arc length along the lifting surface in the 
planes y = constant. 


1 . 5 Method of Solution 

In Ref, 1, it is shown that the distribution of the per- 
turbation aerodynamic potential around a body of arbitrary 
shape is given by the following integral expression: 


4irE4) = 


f 


3(f) 

3n 


1 ^3 Am 

r ^3n^r^^ 


dE 


(1.24) 


where 


E = 0 
E = 1/2 
E = 1 


inside E 

o 

on E 

o 

outside E 

o 



is a surface surrounding the body and . its wake , and n 
is the normal to the surface. 

If the distance between the upper and lower sides of 
the body surface goes to zero (.zero-thickness body) , one 
obtains a lifting surface formulation; 


<j) = 



(-) dl 
r 


(1.25) 


where E extends over the lifting surface and its wake. 


D 


u 


1 ^ 

4it 


(1.26) 


and the subscripts u and 5- stand for upper and lower surfaces, 
respectively. Equation (1.25) shows that the potential can 
be represented by a doublet distribution on the body and on 
the wake. The value of D is constant along streamlines of 
the wake and equal to the value at the trailing edge of the 
wing (Eq. 1.19) . 

The boundary condition, Eq. (1-.8), must be satisfied. 
Using Eq. (1.25), the following integral equation results; 


atj) 



[D (i) ] 


(1.27) 


where 3<j)/9n^ (the subscript zero denotes the control point) 
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is known and given by E. (1.8). 

The sTirface of the wake is assumed to be known (say 
from independent calculations) . 

The numerical solution of Eq. (1.27) will be analyzed 
in detail in this work . 

After Eq. (1.27) has been solved for D, the velocity 
at any point, P, in the field, may be obtained as: 






(-) dE 
r 


(1.28)' 


From Vp, the pressure, as well as a new geometry for the wake 
is obtained. 


1.6 Outline of the Work 

In Ref. 1, the numerical formulation for the integral 
equation describing the distribution of the perturbation aero- 
dynamic potential over a lifting surface has been obtained. 
Expressions for the velocity vector, V, at any point in the 
field have also been obtained. In Section II of this work, a 
summary of Ref. 1 is presented. A description of the itera- 
tion scheme used for obtaining a rolled-up wake geometry, as 
well as the calculation of the nonlinear pressure coefficient 
are added. 

Section III presents results obtained with the lifting 
surface formulation of Ref. 1, shown in comparison with 
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exi'sting theoretical and experimental results. 

The convergence of the solution is illustrated in 
Appendix A. In Appendix B, the convergence of the itera- 
tion scheme is presented. A flow chart and list of the 
computer program implementing the theoretical formulation 
is contained in Appendix C. 



SECTION II 


NUMERICAL POPULATION 


2 . 1 Introduction ‘ 

ThiS' section presents the numerical formulation used 
here, including the wake roll-up iteration procedure and 
the calculation of the pressure coefficient, using the 
linearized Bernoulli Equation, as well as the nonlinear ized 
one. This formulation is an extension of the one of Ref. 1,. 
where wake roll-up is not included. For completeness, the 
formulation of Ref, 1 is summarized here. 

As mentioned in the previous section, the finite-ele- 
ment formulation yields the distribution of the doublet 
strength at the centroids of the lifting surface elements. 

Once this is knov/n, the velocity at any point in the field, 
in particular at the corner points of the wake elements may 
be obtained. These may be used to obtain the geometry of 
the wake . 

In Subsection 2.2, the gradient of Eq. (1.25) is expressed 
in terms of the values of D at the centroids of the elements; 
the boundary condition, Eq. (l.S) is satisfied at the centroids 
of the elements (control points) . In Subsection 2.3, a new 
type of surface element, the hyper boloidal quadrilateral ele- 
ment, first introduced in Ref. 16., is briefly presented, to- 
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gether with the vector expressions for the velocity induced 
by an element at a control -po'int . In Subsection 2.4, the 
iteration scheme used for obtaining the rolled-up wake pat- 
tern is presented. The element grid used .for performing the 
numerical calculations is described in Subsection 2.5. In 
Subsection 2.6, the finite~dif f erence procedure for calcula- 
ting the pressure coefficient in terms of the planform 
geometry is indicated.. 


2.2 Discretization 

The lifting surface and its wake are divided into small 
surface elements. See Pig. 2. Assume that the value of D 
is constant within each element, say it is equal to D (un- 
known) at the centroid of the element a. . Then Eq., (1.27) 
reduces to : 


8c[) 

3 % 



a" 

9n3n 

o 





( 2 . 1 ) 


where N is the number of surface elements on the wing and L 
is the number of elements on the wake. Note that D is con- 
stant along streamlines of the wake and equal to its value 
at the trailing edge or approximately equal to D at the 
centroids of the wing elements in contact with the trailing 
edge. If we impose that the boundary condition, Eq. (1.8) 



- 17 - 


is satisfied at the centroids of the wing surface 

elements the following system of linear algebraic 

equations is obtained: 



( 2 . 2 ) 


( 2 . 3 ) 


( 2 . 4 ) 


( 2 . 5 ) 


( 216 ) 
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••A, , = n. 
hk h 


V, 


hk 


C2.7) 


2,3 Hyperboloidal Quadrilateral Element 

In order to evaluate Egs. (2.3) and (2.6), a typical 
quadrilateral surface element is approximated by a portion 
of a hyperboloidal paraboloid passing through the four cor- 
ner points. This type of surface element is called the hyper- 
boloidal quadrilateral element, introduced in Ref. 16 and 
briefly described here. 

The geometry of a surface element is described in vector 
form as : 


-> ^ , 1 2 

P = P (2.8,) 

1 2 

where 5 ^iid ^ are the generalized curvilinear coordinates 
on the surface elements with the base vectors 









d p 




(2.9) 


The unit normal to the surface is ■ 
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n 



X a. 


X a. 


( 2 . 10 ) 


The surface element is (see Fig. 3) 


da = ja_d^^ x a^d^^ 


a^^ X a2!d5 d^ 


( 2 . 11 ) 


The h’yperboloidal element approximating the real surface 
element is described by the expression (see Fig. 4) : 


P = [1, 71 , Co] 


p. 


p. 


( 2 . 12 ) 


P. 


with 


-1 < C 1 1 
-1 < n < 1 


(2.13) 


where P^ represents the centroid of the element a^^. The coor- 
dinates of the corners of the element are related to P^^, P^, 

?2 and P^ as 
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where (see Fig. 5) : 


-> 





(2.17) 


2.4 Iteration Scheme for Wake Roll-up 

As mentioned in Section I, the wake is initially assumed 
to consist of straight vortex lines- starting at the trailing 
edge of the wing. It was also found that these vortex lines 
should be tangent to the velocity vector and this provides 
the condition for obtaining the rolled-up wake geometry. The 
following iteration scheme is used for aligning the initially 
straight-wake streamlines with the velocity vector: compute 
the doublet strength' distribution at the centroids of the 
elements,., with the wake influencing only the terms 

of the .elements, in contact v?ith the trailing edge. Then 
compute the velocities in the x, y and z directions on the 
wake,' at the corners of the surface elements. Align' segments 
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of .the wake streamlines with- the velocity vector evalua- 
ted at .the upstream segment extremity. (See for example 
Fig. 6, where the position of the point is changed 

according to the velocity at the point on a typical wake 


surface element) 
as follows ; 


The position of the point is changed 


-> -> 

P = P + AP 
pm mm. 


(2.18) 


where 


AP = V jA?l / |v| (2.19) 

and I AP j is the original distance between the points P^^ and 
^mm’ doublet strength distribution is calculated again 

(notice a very small change, due to the new wake geometry) , 
then the wake velocities and geometry are reevaluated. The 
process repeats itself until the difference between suc- 
cessive wake geometries becomes sufficiently small, thus 
indicating the convergence of the scheme (or the fact - that 
the streamlines are indeed tangent to the velocity vector) . 

The iteration scheme described here is hot the best 
possible one. A number of improvements are suggested in 
Appendices A and B. 
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2,5 Element Grid 

The pressure coefficient for the wing is computed by 
using the f inite~dif ference method. In order to properly 
illustrate the scheme, a description of the element grid 
is in order. 

Let c (y) be the chord and b that span of the wing, x 
and y the Cartesian coordinates for the wing at zero degrees 
angle of attack (see Fig. 2). Let 


Y 


c (y) 


e 



( 2 . 20 ) 


Then the parametric form of the wing planf orm equation is : 


^ ’'l.e 



z = 0 for a flat lifting surface (2.21) 


If, in addition, the wing is at an angle of attack, a, 
different from zero, the geometry may be rewritten as: 


X = X cos a 


y = y 
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2 = -X sin a 


( 2 . 22 ) 


Since the potential (doublet strength) varies faster 
near the leading edges and tips of the v/ing, it was found 
convenient to use smaller boxes in these regions and larger 
ones elsewhere. See also Ref. 15- This is accomplished 
by the following transformation: 


,2 

Y = ij; 

g = 1 - (1-0)^ (2.23) 


The boxes have constant si 2 'es in the plane ijj and 0, given 
by: 


ijj = 1/NX 

0 = 1/NY (2.24) 

where NX and NY are the numbers of boxes along the x direction 
and along the semispan, respectively. 


2.6 Pressure Coefficient 

As shown in Subsection 1.4, the linearized pressure coeffi- 
cient ACp is given by 



1^ (A(})) + 0(e) 


(2.25) 
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where s is the arclength on the wing on the plane's y = 
constant. As mentioned before, by solving Eq. (2.2), the 
potential distribution is obtained at the centroids of the 
wing elements. By interpolation, a continuous distribution 
can be obtained. 

The wings used here for the numerical examples are all 
rectangular flat surfaces, for which s = x. The derivative 
of the potential in Eq. (2.25) can be written as; 


3 (A(i)) _ a (A(|)) 

ax " 


aif) 

ay 


ax 


(2.26). 


At any point x^, on an element borderline along the semispan, 
the derivative of the potential, by finite - differences, is: 


a A(j) ^ 1 ^‘^i+1/2 ~ ^‘^i-1/2 


a X . 

1 


^i+l/2 


'^i-1/2 


2/lpT 


(2.27) 


where i ± 1/2 represents adjacent element centroids. on planes 
y = constant. 

The non-linearized pressure coefficient is given by Eq. 
(1.16), reproduced here: 

Ac ■= - 2(1 + v<j> ) • V{A(f))= 2v • V(A^) (2.28) 

P 9- 


Denote by i, j and k the unit vectors along the x, y, z coor- 
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dlnates-and by i , 3 , k the unit vectors along the x, 

w . w w 

y, z coordinates. One can express the velocity V in terms 
of the wing coordinates and in terms of the x, y, y coordi- 
nates as: 


V 


V, + V. : + V_ S 

1 w -2 -"w 3 w 


= Vi + V3+V^ 
X y-^ z 


(2:29) 


and 


V(A(j>) 


9 (A(j)) ^ 
3 X 


3 (Ad)) f 

w ^ — -'w 

3 y 


(2.30) 


Therefore, combining Eqs. (2.28), (2.29) and (2.30), the 
pressure coefficient becomes: 


Ac =2v, iiAil +2v 

P lax 23- 


(2.31) 


On the plane y = 0, a simpler expression may be obtained, 
since, for the symmetric cases considered here 


8 (A4>) ^ 
3 y 


(2'.32) 


Therefore, on y = 0 
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A c =2V- 

P ^ 3 X 


(2.33) 


where is given by 


V, = V • i 
1 w 


V cos a 

X 


V sin a 
z 


(2.34) 


with 


V = 1 + 


Z 3z 


(2.35) 


- 5 —^ are obtained from Eq. (2,5) 

o Z 


Finally, 


^T-i , ^*^a •. ls(A(i)) 

(1 + -T~)cos a - sin a y ~- 

' J 8 X 


(2.360 


where - is computed according to Eqs. (2.26) and (2.27) 

3 X 



SECTION III 


NUMERICAL RESULTS 


3 . 1 Introduc-bion 

As mentioned in the beginning of Sections I and II, 
this v/ork is an extension of ,Ref. 1. The zeroth order for- 
mulation described in Section II was implemented into a com- 
puter code, ILSAWR (acronym from Incompressible Lifting 
Surface Aerodynamics with Wake Roll-up) . ILSAWR performs 
the iteration routine described in Subsection 2.4. The 
v;ay the program is set up, the wake geometry is automatically 
generated, with each rov? of elements along the x-dire‘ction 
having equal lengths. The. Kutta condition is satisfied by 
imposing that the first row of wake elements is tangent to 
the wing. The iteration scheme is performed for the rest 
of the rows only. 

All numerical results presented here were obtained for 
rectangular planar lifting surfaces and all the. graphs show 
results only for the semispan of the wings. 

3.2 Parametric Analysis of the Effect of the Angle of Attack 

A parametric analysis of the effect of the angle of at- 
tack on the wake roll-up is presented here. The case consi- 
dered is a rectangular wing of aspect ratio AR = 8. This 
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value was chosen because of existing results of Ref. 11 
(see Section 3.3). Results. are presented for three values 
of the angle of attack; a = 5°, 10° and 15°. The case 
a = 5° is presented in detail. In Figures 7a, b, c, and d, 
the converged wake pattern for a rectangular planform 
of aspect ratio AR = 8 at an angle of attack a =5°, with 
an element grid having NX == 4, NY = 10, with the length 
of the wake elements Ax^ = .5c is plotted in great detail 
for 10 chord lengths behind the trailing edge. Figures 
7a and 7b show the rolled-up wake plotted at stations 1 
through 10 chord lengths behind the trailing edge. Figure 
7c is a side viev/ of the rolled-up wake (the vertical scale 
is enlarged) , showing the vertical displacement of the stream- 
lines. The numeration system for the streamlines is also 
shown, with streamline number 1 being at y = 0 and the last 
streamline starting at the wingtip. Figure 7d shows> a top 
view of the rolled-up wake behind the wing, with the side 
displacement of the streamlines visible. The streamline 
numeration system is clearly shown here. 

It may be noted that the analysis of convergence (pre- 
sented in Appendix A) indicates that the solution is close 
to convergence, although improvements appear to be desirable 
at the trailing edge, especially near the wingtip. 

Results for a = 5°, 10° and 15° are presented in Figures 

8 and 9 , for a rectangular planform of aspect ratio AR = 8 , 

with an element grid of NX = 4, NY = 10, Ax = ,5c. Figure 8' 

w 
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shows the effect of the angle of attack on the wake roll- 
up, plotted at 5 chord lengths behind the trailing edge. 

The wake displacement becomes more pronounced as a increases 
from 5® to 10°' and 15°. The effect of the angle of attack 
is shown also in Figure 9a, b and c, where streamlines 
1, 10 and 11 are plotted in a side view. 

The analysis of the convergence of the iteration scheme 
is presented in Appendix B. 

3.3 Comparisons with Existing Results 

In order to assess the validity of the method, a number 
of comparisons with existing results are presented here. 

3.3.1 Comparison with the Artificial Viscosity Method of 
Bloom and Jen 

Figure 10 presents the wake roll-up for a rectangular 
planform of aspect ratio AR = 8 at an angle of attack a = 6.25°, 
for an element grid of NX = 4 , NY = 10 , with the wake elements 
length of Ax^ = . 5c . Converged wake patterns are shown at 
stations 1, 5 and 9 chord lengths behind the trailing edge 
and the results of the present method are compared with the 
artificial viscosity results of Ref. 11. In Ref. 11, the 
lift coefficient was C^. = 1 and no angle of attack was sped- 
fied. Therefore, the lift coefficient per unit angle of 
attack, , was evaluated with -the present method and the 
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angle of attack' was found according to 



The lift coefficient per unit angle of attack was found to 

be = 9.17'4. For this value of the C,. , the value of the 

La. La 

angle of attack which gives a lift coefficient of 1 is a = 6.25°. 

3.3.2 Comparison with the Experimental Results of Chigier 
and Corsiglia 

In Ref. 8, the position of the vortex centerline is 

experimentally determined as the locations where the tangential 

velocity is zero. The results Of Ref. 8 (Chigier and Corsiglia) 

have been obtained for a rectangular wing of aspect ratio 

AR = 6 , at an angle of attack of a = 12°. For the present 

method, there is (as yet) no precise way for determining the 

location of the vortex centerline. The last streamline is 

taken to represent the vortex centerline for the planform with 

an element grid of NX = 4 , NY = 10 and Ax = .5c. Figure 11 

w 

results obtained with the present method, compared with the 
ones of Chigier and Corsiglia. 

3.3.3 Comparison with Results of Shollenberger 

As mentioned in Section I, Ref. 12 (Shollenberger) uses 
a three-dimensional potential method and an iteration proce- 
dure to obtain the rolled-up wake. The wing planform used 
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has an aspect ratio AR = 6 and it is at an angle of attack 
a = 10°. The results obtained with the present method, in 
comparison with the ones of Ref. 12, are shown in Figure 12. 
The wake geometries are plotted for 1, 2, 3 and 4 chord 
lengths behind the trailing edge. 

I 

3.4 Pressure Coefficient 

In Subsection 2:6, the finite-difference procedure used 
in calculating the pressure coefficient was described in de- 
tail- The results obtained by using the linearized and non- 
linear Bernoulli Equations with and without wake roll-up are 
presented here. Table I shows the values of Ac^ at y = 0, 
linear and nonlinear, with straight and rolled-up wakes. 

The results are obtained for a rectangular wing with aspect 
ratio AR = 8, at an angle of attack a = 5°, with an element 

grid of NX = 7 , NY = 7 , with Ax = . 5c . Figure 13 shows a 

w 

plot of the pressure coefficient presented in Table I. Note 
the negligible effect of the wake roll-up on Ac^. However, 
as previously mentioned, the wake roll-up is believed to have 
an important effect in the case of wing-tail interaction. 

Finally, Figure 14 presents the potential distribution 
at the trailing edge of the wing, Ad),_ _ , for the same plan- 
form as the one used in Figure 13 . It can be seen from Figure 
14 that at y = 0, 8(A^)/3y = 0. The effect of the wake roll- 
up is negligible. 
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x/c 

Linearized 

Ac ■ 

P 1 

Nonlinear 

Ac 

P 


Straight 

Wake 

Rolled-up 

Wake 

Straight 

Wake 

Rolled-up 

Wake 

.055 

.8800 

.8706 

, .8699 

.8600 

.136 

.4523 

.4530 

.4^471 

.44-70 

.258. 

.2799 

.2860 

.2767 

.2827 

.421 

.1921 

.1921 

.1899 

.1899 

.624 

.1276 

.1274 

.1261 

.1259 

.868 

.0732 

.0729 

.0723 

.0720 


Table 1. Pressure Coefficient at y = 0, for a rectangular 
wing planf orm of aspect ratio AR = 8 , at angle of attack 
a = 5®/ with element grid of NX = 7, NY = 7 and Ax = . 


5c 




CONCLUDING REI^RKS 


A method for analyzing the wake roll-up has been des- 
cribed and numerical results have been presented. At this 
point, it might be interesting to quote Ashley and Rodden 
(Ref. 17) from their review on wing-body aerodynamic -inter- 
^action: "It should be evident from the foregoing all too brief 
account of interaction theory that it is both a complicated 
subject and one in which computer automation is more nearly 
in a state of revolution than of evolution. Within a few 
years, programs should be available that will solve the 
linear potential equation ^ with boundary conditions satis- 
fied by placing appropriate discrete singularity elements at 
a close approximation to ail the true wing and body surfaces. 
The following 'nonlinearities' will be i.ncluded: pressure 
velocity relations such as" the nonlinear Bernoulli Theorem; 

"boundary conditions that partially account for x - velocity 

* 

perturbations . . . ; v/akes trailing streamwise from the actual 
positions of trailing edges; and/or estimates of self-deforma- 
tion of wing wakes as they affect aft tail surfaces and the 
like. " 

All the nonlinearities mentioned above (with the exception 
of the wing- tail interaction and the zero- thickness limitations 
of lifting-surface theory) have been included in the present 
work.' The only approximations introduced are numerical" ones, 
and they are negligible, as the convergence analysis 
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indicates . 

Finally, the main innovations and advantages of the mehod 
are discussed. First, the method is based upon an exact 
(rather = than discrete) formulation. Only numerical approxi- 
mations are introduced (other methods use approximate physical 
models such as discrete vortices) : this implies that the 

formulation is apt to refinements (first-order finite-element 
representation for D is now under • investigation) . Second, 
the wake is represented as doublet distribution: this implies 
that the method may be extended to steady and unsteady, sub- 
sonic and supersonic- aerodynamics around complex configura- 
tions, in a relatively straightforward method, using the for- 
mulation of Ref. 18. Third, the convergence, of the solution 
is exceptionally fast (as is the more general method of Ref. 

18). Fourth, the method is relatively fast; the results for 
NX = 4, NY = 7, = 10 require 3 minutes of C.P.U. time 

per iteration on the I.B.M. 370/145 computer of Boston 
University. Finally, the convergence of the iteration scheme 
is already good, although considerable improvements can be 
obtained by using alternative, more sophisticated iteration 
schemes which are now under investigation. 

Most of the theoretical results on wake roll-up are of a 
rather recent origin (from 1973 onward) and comparisons with 
experimental results show that some refinements of the mathe- 
matical model are still in order. Viscosity effects, thickness 
effects, aerodynamic interaction still remain to be accounted for. 
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Figure 2 . 


Lifting Surface and Wake Geometries. 
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Figure 4. Hyper boloidal Quadrilateral Element. 
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Figure 7a. Converged Walce Pattern, gor a I^ecta,agular Wing pianfprm 

of AR = 8, with a = 5°, Element Grid with NX = 4, NY = 10, 
Length of Wake Elements Ax^ = .5c, Plotted for 10 Chord 
Lengths Behind the Trailing Edge. Continued on Next Page. 






LIFTING SURFACE 




Figure 8, The Effect of the angle of attack, a, on the Wake 
Roll“up for a Rectangular Lifting Surface of 
AR = 8, Plotted at 5 Chord Lengths Behind the 
Trailing Edge. The Element Grid has NX = 4 , NY = 
and Ax = ,5c. 
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Figure 10. Wake Roll-up for a Rectangular Lifting Surface of AR = 8, 

Angle of Attack a = 6.25°, Element ' Grid with NX = 4, NY = 10, 
Ax = , 5c and Comparison with the Results of Ref. 11. 






Figure 12, Wake Roll'^up for a Rectangular Lifting Surface of AR = 8, 
at a = 5°, with NX = 4, NY = 10, Length of Wake Elements 
= ,5c and Comparison with Results of Ref, 12. 

Continued on Next Page, 
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APPENDIX A 


CONVERGENCE OP SOLUTION 


In -this Appendix, a numerical study is performed on the 

influence that the parameters NW (the number of wake elements 

along the x direction ) , Ax and NY have on the convergence 

w 

of the solution. The only case presented here is relative to 

a rectangular wing planform of aspect ratio AR = 8 at a = 5°. 

In Figures Ala, b and c, the effect of the length of 

the wake elements on the wake roll-up is shovm. The element 

grid for the planform has NX = 4 and NY = 10 . The length 

of the wake elements Ax is allowed to vary from .5c to .75c 

w 

and Ic. Rolled-up wake patterns are plotted for stations 
located at 3, 6 and 9 chord lengths behind the trailing edge. 
Note that, as Ax increases, the wake pattern becomes "larger", 
as it is easy to see from Figure Al. Note also that the 
difference between wake patterns for various values of Ax^ 
becomes smaller as the distance from the trailing edge in- 
creases. In Figure A2 one might find an explanation to this 
difference, as well as a suggestion for the improvement of 
the numerical model. In this figure, the same rectangular 
planform v/ith an element grid of NX = 4, NY = 10 is used. The 
figure shows streamlines (counted from the line of symmetry of 
the wing, the mid-line included) numbers 1, 10 and 11, plotted 
for values of Ax of .5c, .75c, and Ic, for 10 chord lengths. 
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-The first streamline shov/s remarkable closeness (on this 
enlarged vertical scale) for the various Ax^. The difference 
increases as we approach the wing-tip streamlines. Note 
that the streamlines are approximately parallel; the dif- 
ference between them is due to the fact that, by imposing 
the -Kutta condition, the first row of wake elements lies 
.-in the same plane as the wing, and since A’X^ -varies, the 
streamlines will start at .5c, .75c, and Ic behind the 
trailing edge. Also, the downv;ash is larger in the vicini- 
ty of the trailing edge and decreases as we move farther 
behind. Therefore, the wake slopes can be expected to be 
larger in the vicinity of the trailing edge. Note the 
sharp jump between the first element and the next in 
streamline number 11. It may be worth noting that, since 
the streamline displacemenb is obviously influenced by the 
length of the wake elements, we might obtain a smoother pro- 
file in the vicinity of’ the trailing edge by using smaller 
elements in this region, for one or two chord lengths. This 
remains to be implemented. 

Next, consider the influence of NW. If the number of 
wake elements is increased, it is observed’ that the newly 
added rows of elements have no effect v?hatever on the wake 
roll-up of the previous ones. 

The effect of the number of wake strips on the wake roll- 
up is shown in Figure A3 , for the rectangular wing having 
an element grid of NX = 4, with NY varying between 7 and 10,- 



-60- 


= ,5c and NW - 11. All cases are converged and lie prac- 
tically on the same line. 



WING TRAILING EDGE 


Figure Ala. Influence of the Length of the Wake Elements Ax 
on the Rolled- up Wake fer a Rectangular Wing 

of AR = 8, for a = 5°, at a Station Situated 
3 Chord Lengths Behind the Trailing Edge. The 
Element Grid has NX = 4, NY =' 10. 





Figure Ale. Influence of the Length of the WaKe Elements 

on the Rolled-up . _ Wake for the Wing of Figure A1 
Plotted at a Station Situated 9 Chord Lengths 
Behind the Trailing Edge. 















APPENDIX B 


CONVERGENCE OF ITERATION SCHEME 

In this. Appendix, an analysis of the convergence of the 
iteration scheme is presented, for a rectangular wing of AR = 8, 
at an angle of attack a = 5°, with an element grid having NX = 

4, NY = 10 and Ax = ,5c. Figures Bla, b, c, d show the evo- 
lution of the rolled-up -wake pattern through successive itera- 
tions until convergence is reached. The plots are for stations 

I 

at 1, 2, 5 and 10 chord lengths behind the trailing edge. Fi- 
gures B2a, b, c, d show- the evolution of the wake streamlines 
numbers 1, 9, 10 and 11 through successive iterations until 
convergence, plotted for 10 chord lengths behind the trailing 
edge . 

A common feature of Figs. Bl and B-2 is that the plot& 
of the initial iterations show very large displacem.ents of 
points on the wake. The largest displacement takes place 
near the, wing-tip and far behind tlie trailing edge. Con- 
vergence is attained faster near the trailing edge and the 
rate of convergence decreases as we move from the wing root 
toward the wing-tip . 

The computation time required to obtain the convergence 
of the iteration scheme (described in Subsection 2.4) for a 
wake having 210 elements is of the order of one hour and 20 
minutes on Boston University's IBM 370/145 computer. A num- 
ber of improvements of the present iteration scheme can be 
tried. 

First, since the calculated potential distribution on 
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the wing • is essentially the same with a straight wake as 
well as with a rolled— up one, the potential distribution could 
be computed for the straight wake and then recomputed for 
instance every fifth iteration. This should lead to some 
savings in computational time. Second, a much better itera- 
tion scheme can be used (suggested by the plots of Fig^pres 
•B1 and B2) , This scheme should converge' Ttiuch faster ;t'han 
the one used in this paper and account for significant time 
savings. The first iteration should only change the po- 
sition of the second row of boxes on the wake; (the first 
one is kept tangent to the wing plane according to the way 
the Kutta condition is satisfied) the rest of them will 
have the same y and z coordinates as the second row. Only 
the velocities at the influencing corners are calculated. 

The third row of corners .should be .realigned according to 
the velocities at the second row; the rest of the boxes will 
have the same y and z coordinates as the second row. The 
process should be repeated until convergence is .reached every- 
where . 
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# 10 , 11 , 12 
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Figure Blb^ Evolution of the Rolled-up Wake Pattern 

Through. Successive Iterations., a.t a Station 
Situated 2' Chord Lengths Behind the Trailing 
Edge, for the Wing Planform of Figure Bla 
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F5.gure Bid. Evolution of the Rolled“Up ‘Wake Pattern 

Through Successive Iterations, at a Station 
Situated 10 Chord Lengths Behind the Trailing 
Edge, for the Wing of Figure Bla. 










APPENDIX C 


FLOW -CHART AND LIST OF THE COMPUTER 
PROGRAM ILSAWR 


Cl« Flow Chart of " Computer Program ILSAWR 
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C2, List of ^Computer Program ILSAWR 


I-V G LEVEL 21 MAIN DATE = 75226 05./22/L 

COMMDN/ZZZl/NX ,NY 7 NZ, NW t R E'^L EN , SP AN » KSYMM Y, K SYMMZ 7 N SYMMY , NSYMMZ 
CO»MON/IZl2/T AU, ALFA, T ANGLE , TANGTE 7 CHORD s NTOT AL , UMACH 
COMMON/ ZZZ 8 / A A { 2500) , SGURCF (250) , SINAGC 7 CCS ABC 7 ALF'ASC 
COMMON/ 2 ZZH/VHKX (250 0) 7 VHK Y { 2500 ) , VHKZ ( 2500 ) 7 VKX ( 250) , VK Y{ 250) , 
lVKZ(-250> 

COMMQN/CONTR/NITER 
DIMENSION ITCCNT(IOO) 

DO 10 I=l7l5 
10 ITC0NT{I)=1 
ITCONK 1)=2 
ITCONTl 5)=2 
ITC0NT(10)=2 
ITCONK 15)=2 
ITCONT(20)=2 
ITC0NT(25)=2 
ITCONK 30} =2 
ITC0NK35)=2 
ITC0NT(40)=2 
ITC0NK45)=2 
ITC0NT(50)=2 
CALL INITIA(l) 

CALL PRINTA(5) 

CALL GEOMET 
CALL VEC123 
C CALL PRINTA{3) 

00 1 NITER=l,l2 

I F(NITER.Eq.ll) ITCONK NITER )=2 
IF(NITER.EQ.12) ITCONK N I TER ) = 2 
IFnTC0NT(NlTER),NE,2 )G0 TO lOOO 
CALL COEFF 
C . CALL PRINTB(4) 

T0L=0.001 

CALL GELG{ SOURCE, AA.N-TOTAL,! 7 TOL, I'ER')' 

CALL PRINTBU) 

1000 CONTINUE 

GALL VELMM 
C CALL VELAUX 

CALL ITER 
1 CONTINUE 

STOP 
END 
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SURROUTINE INITIAfK) 

COMMON/ ZZZl /NX, MY, MZ, NW , R FFL FN , SP AN , K SYMMY, KSYMMZ , NSYMMY , N.SYMMZ 
COMMON/ ZZZ2/TAUt ALFA, TANGL E , T ANGT E, CHORD, NTOT AL , IJM'ACH 
C0MM0N/ZZZ3/YK( 3, 11,11,2) 

C0MM0N/ZZZ6/XPC(250 ), YPC (250) , ZPC(250) 

C0MM0N/ZZZ7/XPU250), YPl (250) , ZPl ( 250 ) , X°2 ( 250 ) ,YP2(250) , 

1ZP2(250) , XP3( 25G) , YP3 (250) ,ZP3(250) 

C.OMMON/ Z.2 28/ AA(25 00 } , SOUP CF ( 2 50 ) , S IN A SC , COS ARC , ALF A3C 
C0MMQN/ZZZ9/XPP{250) ,Y?Pi 250) ,ZPP (250) ,XPM(250) ,YPM(250) 

I', ZPM-( 250) ,X'MP (2 50 ), YMP(250) tZMP (250) , XM;M( 250) ,YMM( 250) , 

IZMM (250) , IWAKE(250> 

CnMM0N/ZZZlO/JNXS(25O) ,NXWAKE ,WAKEIN 
■G.0 T0-(1,2,3,4),K 
1' C'ON-TINUE 
NX = 7 
NY=7 
N’Z=1 

NX-WAKE=ll 
WAKEIN=.5 
MW=1 

H MEANS THE GEOMETRY OF THE PROBLEM IS SYMMETRIC 

-i MEANS THE GEOMETRY OF THE PROBLEM IS ANTISYMMETRIC 

0 MEANS THE GEOMETRY QP THE PROBLEM IS NEITHER SYMMETRIC NOR ANTI 

IF KSYMMZ «NE. 0 ,THEN NZ=I (EXCEPT FCR GROUND EFFECT) 

KSYMMY=+-1 
KSYMMZ=0 

NTOTAL=NX*NY*NZ*NW 
I F { KS YMM Y . E Q . 0 ) NT QT AL= NT G T AL * 2 
NSYMMY=1 
NSYMMZ=l 

IF.{KSYMMY.NE.0)^'SYM^'Y=2 
IF(KSYMMZ.NE.0)N$YMMZ=2 
UMACH=0,0 
REFLEN=1. 

TAU^.OO 
SPAN=8. 

ANGLA=0. 

ANGLB=6, 

ALFA^5. 

ALFAR==ALFA’!'3, 14159/ 130. 

SINALF=SIN(ALFAR) 

COSALF=COS(ALFAR) 

alfabc=o. 

ALFRBC=ALFA3C^3. 14159/180. 

SINABC=SIN{ALFRBC) 

COSABC=COS{ALFRBC) 

BETA = SQRT( 1.-UMACH={=UMACH} 

XLEZ=-1. 

XTFZ=Oc 

CHnRD=XTEZ-XLEZ 
XLEZ=XLEZ/(REFLEM>>^BFTA) 

XT'"Z=XTEZ/(REFLEN’i'BFTA) 

SPAN=SPAN/REFLEN 
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HFSPAN=, 5*SPAN 
XLED=ANGLA/BETA 
XTEP = ANGLB/Bt:TA 

TAU8AR=TAU=«".75*SQRT{3..)*{ XTEZ-XLEZ) 

RETURN 

2 CONTINUE 
RETURN 

3 CONTINUE 
nxx-1 . /NX 
DYY=1./NY 
NXP=NXfl 
NYP=NY+1 
DO 33 IX=l,NXP 
DO 33 IY=1,NYP 
DO 33 1Z=1 tNZ 
XX=I IX-1 )*OXX 
YY=nV-l)«nYY 
CSI=XX^XX 

THIS IS FOR A UNIFORM Y-MESH 
ETA^YY 

THIS IS FOR A NCNUNIFORM Y-MESH 
ETA~l,-( I .“YY)^^2 
Y=HF$PAN*ETA 

THIS IS A SEMI-ELLIPTICAL WING PLANFOR'M 
XLE=-CHGRD^’'SQRT (I ,-(Y*Y ) / { HFSPAN*HFSPAN ) ) 
THIS IS A RECTANGULAR WING PLANFORM 
XTE=XTEZ 
XLF=XLEZ 

X0=XLE4-{ XTE-XLE)*CSI 
IF{ IZcEQ.l)SIGNZ=.^-l 
IF( IZ.E0.2)SIGNZ=-l 

ZO=SIGNZ*TAUBAR=^XX« a .-CSI )*SGRT( l.-ETA**2) 
X^fXO’l'COSALF FZO^S INAL F 
Z=-XO*SINALF+ZO«CCSALF 
YK{ Ij IX, lY, IZ) = X 
YK(2,IX,IY,IZI=Y 
YK(3, IX, IY,IZI=Z 
33 CONTINUE 

RETURN 

4 CONTINUE 
RETURN 
END 
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subroutine GECMET 

THIS SUBROUTINE IS FOR- QUADRILATERAL ELEMENTS 

C0MM0N/ZZZ1/NX,NY,NZ, NW, REFLEN , SP AN t KSYMNY, KSYMMZ » NSYMMY , NSYMM7 
CQMNON/ ZZZ2/T A'J , ALFA, TANGLE, TAN GTE, CHORD , NTOTAL , UM ACH 
CCYMON/ZZZ3/YKn, U ,1 1,2 ) 

C0NM0N/ZZZ6/XPC (250),YPC(25D) ,ZPC(250) 

CGMMON/ZZZT/XFl (250 ), YPl (250) ,ZP1 (250), XP2( 250 ) , YP 2 ( 25 O') , 
1ZP2(250),XP3(250),YP3(250),ZP3(250) 

C0MM0N/ZZZ9/XPP (250 ), YPP { 250) ,ZPP (250) , XPM{ 250) ,YPM(250) 
1,ZPM(250),XMP(250),YMP(250), ZMP (2 50 ) , XMi«{ 250 ) , YNM( 250} , 

IZMM (250) , I WAKE { 250) 

COMMON/ Z ZZ 10 / JNXB( 250) , NXWAK F , WAK E IN 

INDEX(JW,JX,JY,JZ,MW, MWX,MWXY)=JW<-MW*(JX-1> fMWX*(J Y~1 )' i--MWXY*'( JZ —1 ) 

NWX=NW*NX 
NWXY=NWX*NY 
NWXYZ---NWXY*NZ 

CALL INITIAOJ 

DO 200 IX=1,NX 
DO 200 IY=L,MY 
00 200 IZ=1,NZ 

+ + 

IW=l 

IN0=INDEX( IW, IX, lY, 1Z,NW,NWX,NWXY) 

IF( IZ.E0.2)G0 TO 906 
IXMM=IX 
lXPM-IXi-1 
IXPP= IX + 1 
IXMP=IX 
I YMM=IY 
I YPM= lY 
1YPP= lYM 
IYMP=I Y^-l 
IZMM= IZ 
J ZPM=IZ 
IZPP=IZ 
IZMP= IZ 

+- 

GO TO 999 
906 CONTINUE 
I XMM=TX 
IXMP=IX 
IXPP=IX^l 
IXPM=TXH 
IYMM=IY4-1 

iymp=iy 

IYPP=IY 
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r-YP-M=IY+l 

I,ZMM=I2 

IZMP^IZ 

IZPP=IZ 

IZPM=IZ 

999 CONTINUE 
C 

XPPnND)=YK(l,TXPP,TYPPtI2PP) ' 

YPP{IN0)=YK(2»IXPP»IYPF, IZPP) 

2PD UMH) = YK( 3» IXPP, lYPP, IZPP) 

XPM(TNO)=YK( 1 , IXPY, lYFYt TZPM) 

YPM{ IND)=YK(2, IX^P, lYPM, I ZPM) 

ZPM(IN:D)=YK(3,IXPP, lYPP, TZPW) 

XMP (IND) =YKU , IXMP, lYYP, I ZMP) 

YMP{ IMH) =YK{ 2f IXPP, lYMP, TZMP) 

ZMP{IND)^YK(3, IX'-^P, lYMP, IZMP) 

X.MMtINO)-YKU, I XMK, I VMM , I Z^1M ) 

YMMn.MD)=YK{2,IXMP, lYMMt IZMM) 

ZMM{ JNO)=YK(3t I Xwp, IYMM» IZMM) 

C WR1TEC6, 199) IND,XPP(IND) , YPP {TND )■, ZPP ( INOl t XPM ( INO) tYP^H I Mn ) 

C I, ZPM{ IND) , XM?‘( IKO) T YMP( IMD) , ZMP UMD) , X^M( TMD) ,YMM ( INH) , ZMM( IND) 

199 FORMAT!/ ’ IND=’ , 1 2 , / • PP » , 3 X, 3P 10. 4/ ' PM * , 3X ,3 F 1 0. 4/ ‘ MP » , 
13Xt3F10.4/'MM' ,3X,3F10.4) 

lWAKFnMr)} = 0 

IF{ IX.EQ.NX) !WAKE(IND)=L 

200 CONTINUE 
IFCKSY'^MY.NE.OIGO TC 701 
DO 300 IR=ItNWXYZ 
IL=TRi-NWXYZ 
XPP{IL)=tXMP( IR) 

•XMP{ IL) = fXPPUR) 

XPM(IL)=+XMM( IR) 

XMM( IL) = -»-XPM( IR) 

YPP{ IL)=-YMP{ IR) 

YMP( IL)=-YPP{ IR) 

YPM( IU=-YMM( IR ) 

YMM(IL)=-YPM( IR) 

ZPPnU= + ZMP( IR) 

ZMP( IL) = +ZPP{ !R) 

ZPM(IL)=+ZMM( IR) 

ZMM( IL)=+ZPM(TR) 

TW4KE( IL )=IWAKE ( IR) 

300 CONTINUE 
701 CONTINUE 
C 

C -POINTS FOR WAKE 
C 

NWTOT=NXWAKE*NY 
no 1 IY=1 tNY 
DO 1 IX=1,NXWAKE 
JNXW=IXf (IY-1 )«NXWAKE 
JNXB(UNXW)-=IY’S=NX 
1 CONTINUE 

00 10 IY=1,NY 
. DO 10 IX=1,NXWAKE 

1 1= IY*NX 

IND'=NTOTAL+IX+( 1Y-1)*NXWAKE 
FACT0R=1. 


EEPRODUCIBILrrY OP THE 
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IX.f=Q.NXWAKP)FACTnR=100. 

XMM( INO) = XPM{ II ) + WAKFIN*( iX-n 
XMP{ IND)=XPP( II )+WAKE-IN*nx-l) 
XPP( I.ND) =XPP ( II ) 4-W'AKE IN*IX*FvA.CTOR 
XP.M{ IND) = XPM.( II )+WAKEIN*lX«FACTOR 
YM-M (IND) =YPM { 1 1 ) 

Z«M( IND)=ZPM( IT ) 

YMP( INO)=YPP( II ) 

ZMPdNni -=ZPP( II ) 

YPP{ IND)=YPP('I I ) 

ZPPnND)=ZPP(II} 

YPN{ INT3 =YPM( II ) 

ZPM( INO)=ZPM{ II ) 

10 CONTINUE 
RETURN 
END 
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SUBRDUTINE VFC123 

COM-MDN/ZZZl/NX'tNY ,NZ, NW , REFLEN , SP AN, KSV , KS YMMZ , NSYMMY» NSYMMZ 
C0M.M0N/ZZZ2/TAU, ALFA, tangle, TANGTE, chord, NTOTAL fUM AC H 
CGMYON/ZZZ3/YK{3, 11 ,11,2) 

C0MMnN/ZZZ6/XPC{2 50), YPC (250) ,ZPC(250')' 

COMMON/ZZZ7/XP1 (250) , YPl (250) ,ZP1(250)-, XP2(250) ,YP2(250) , 
1ZP2(250),XP3(25C),YP3(250),ZP3(250) 

COMMON/ ZZZ9/XPP (250),YPP(250),ZPP(250),XPM(250),YPM(250) 

1, ZPM(250) ,X.MP(250), YMP(250) , ZMP (250) , XYM( 250) ,YMM{250) , 
1Z.MM{250) ,IWAKE(250) 

DO 200 IND=1,NTGTAL 

XPC( rNO) = (XPP (IND)^-XPK{ INC) ^XMp (1 ND ) +XNM { IND) )/4. 

YPC( IND) = { YPP UNO )■^YPM( I NO fYMP ( IND) -f-YMM( IND) )/4, 

ZPC( IND)=(ZPP{TNn)fZPM( INO+ZMP ( IND) +ZM.M( IND) )/4. 

XP1{IND) = (XPP(I^'D) + XP^'( INO-XMp(IND)-XMM( IND) )/4, 
YP1(IND)={YDP{IND)+-YPM{ IND) - YMP ( I NO )- YMM ( IND ) ) / 4 . 

ZP1(IND) = (ZPP(IN0) + ZPM( INO-ZMpnND)-'ZMM{ IND) )/4. 

XP2{IN0) = (XPP(IND)-XPP( INC) i-X PP ( I NO ) -XMR ( IND) )/4. 

YP2( IND*)= (YPP (INDJ-YPM( INO) J-YMP ( I ND ) - YN.M ( IND ) )/4. 

ZP2{IND) = ( ZPP {IND)-ZPM( IND) 4-Z MP ( I NO )-ZYM ( IND) )/4. 

XP3(INn) = (>:PP!IND)-XPH( I ND ) -X PP ( I NO ) f XN5M ( IND) )/4. 

YP3( IN0) = (YPP (TN0)-YPP( INC)-YMP( INO) +YMM( IND) )/4. 

ZP3( INn) = (ZPP(IND)-ZPY( INO)-ZPP(IND)i-ZMM( IND) )/4, 

200 CONTINUE 
RETURN 
END 
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SU5ROUTTN6. PR INTA ( KPR I NT ) 

COMMON/ ZZZI/NXt MY ^NZ, NW , RE^L CN , SPAN , KSYMMY, KSY-MMZ , NS YMM Y» NSYMMZ 
COMMON/ZZZ2/TAU » ALFA, TANGLE, TANGTE, CHGRO,NTOTAL, UM AC H 
COMMON/ ZZZ3/YK( 3, 11,11,23 
COMMON/ Z Z Z6 /XPC { 250 ),Y PC (2503 ,2PC1250) 

COMMON/ZZZ7/XP1(250),YP1(250) 5 ZPl { 2 50 ) , XP2 { 250 ) ,Y‘P2(250) , 

IZP2 (2503 ,XP3( 250) ,YP3 (250) , ZP3{250) 

eOMMnN/ZZZ8/AA(2500) , SOURCE ( 2 50 3 , S I.NABC , COS ABC , AL F ABC 
C0MM0N/ZZZP/XPP(250 3 , YPP ( 25 0 3 , ZPP ( 2 50 ) , XP M { 2503 , YP V { 2 50 3 
1, ZPM(250) ,XMP(250 3 ,YMP(250) , ZMP i 250 3 , XMM ( 250 ), YMM ( 250) , 

IZMM( 2503 , IWAKE( 2505 
COMMON/ ZZZ 10 /JNXB (250) ,NXKAKE, WAKE IN 
COMMON/CONTR/NITER 
NTP=NTOTALU 
NTBW=NTOTALi-NY*MXWAKE 
MY4=4*( NY-1) 

GO T0(1,2,3,A,5,6),KPRINT 

1 CONTINUE 
RETURN 

2 CONTINUE 
NXP=NX+1 
NYP=NY+I 

DO 35 1Z=1tNZ 
00 35 IY=1,NYP 
DO 35 IX=1,NXP 
DO 35 J=l,3 

WRITF (6,2 500) J, IX,IY, IZ,YK( J, IX, lY, IZ) 

2 500 FORMAT! ’YK( * , I 1 , • , * , II , ' , ’ , 1 1 , • , * , 1 1 , ' ) = * ,Fi5.6) 

35 CONTINUE 

RETURN 

3 CONTINUE 
WRITF(6,400 ) 

400 F0RMAT(2X, • INC ,4X, «XPC* ,7X, ‘ YPC’ ,7X, * ZPC ’ , 7X , » X P 1 ’ , TX • YPl ' , 7X, 
1’ ZPl* ,7X, »XP2»,7X, ‘ YP2’ , 7X, • ZP2» ,7X, 'XP3» ,7X,* YP3» ,7X*ZP3« ) 

DO 45 I=1,NT0TAL 

45 ■ WRITE (6,500) I ,XPC ( I ) , YPC ( I 3 , ZPC { I ), XPU I) ,vp 1 n ) , ZPK I ) , XP2( I ), 

1YP2( n, ZP2( I ) ,XP3( I ), YP3( n ,ZP3(I ) 

500 FORMAT! IX, I3r 12F10-5) 

RETURN 

4 CONTINUE 
RETURN 

5 CONTINUE 
WRITE(6, 550.) 

550 PQRMAT(//2X, » SPECTPICATIONS OF THE PROBLEM*/) 

WRITE(6,555)NX,NY,NZ, NW,NTOTAL , KSYMMy, KSYM mZ, RE^LEN, SPAN, TAU, 
lALFA, ALF ABC, UMACH,NXWAKE, WAKE IN 

555 FORMAT ( 2 X, 'NX=* ,I2/2X,*NY=’ , I 2/2X , * NZ= • , I2/2X,*NW=’, 12/ 

X2X, ‘NTOTAL = * , I3//2X,,* KSYMMY= • , I 2/2X, * KSY-MMZ=* ,12// 

12X, ’ reference L ENGTH= » , F6.2/2X, »SPAN/REF LENGTH =*,F6.2/ 

12 X, ' THICKNESS^* ,F9.5//2X , ' ALFA=* , F7.3/2X, ' ALFABC=* ,F7.3// 
12X,*MACH NUMBER = * , F7 . 3/ / 2X , * NXWAKE= • , I 3 , // 

I2X , * WAKEIN-=’ ,F7,3) 

WRITE (6, 556) TANGLE, TANCTE , CHORD 

556 FORMAT (2X, * TANGL E= ' , P6 , 2/ 2X , * TANGTE= • , F6 , 2/ /2X , *CHORO=* ,F6.2//) 
RETURN 

6 continue 
RETURN 
END 
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200 
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770 


771 
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775 

776 
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77 


SUBROUTINF P R TNT3 ( KPR T NT J 

COMMON /ZZZl /NX t NY »NZ, NW t REFLFN , SP AM , KSYMMY» K SYMM-Z » NSYMMY, NSYMMZ 
C0MM0N/ZZZ2/TAU , ALFA, TANG l.F, TAMGTF, CHORO , NTOTAL , UMACH 
COMMOM/ZZZ6/XPC{2 50), YPC( 250) ,ZPC(250) 

C0MMnN/Z,ZZ7/XPl ( 250) , YP1(250),ZP1(250),XP2(250),YP2{250), 
1ZP2(250) ,XP3 ( 250) ,YP3 (250 ) , ZP3 (250) 

C0MM0N/ZZZ3/AA( 2500) , SOURCE (250) , S INA8C , COS ABC, AL FABC 

N'WX=NWX^NX 

NWXY=NWX'=^=NY 

NWXYZ=NWXY=!^NZ 

NY4=4*(NY-1 ) 

GO 70(1,2,3,4,5,6,7), KPRTNT 

CONTINUE 

WRITE(6,100) 

FORMAT( ///2X, ’THE DISTRIBUTION OF THE DOUBLET STRENGTH DH») 

INDFIM^O 

I PRINT=0 

DO 25 ISYMY=1,2 

IF( ISYMY.EQ.2.AND.K$YMMY*NF.>0)GD to 25 
!’={ ISYMY.EQ,! ) WRITE (6,120) 

F0RMAT(//5X, ’RTGHTHAMC SIDE’) 

I F { ISYMY . EQ« 2 ) WR ITE ('6 , 1401 
F0RMAT(//5X, ' LEFTHANO SIDE') 

00 25 IZ=1,NZ 
I^^DFIN=IN^FIN^NWXY 
IPRTMT=IPRINT+1 
INO=NWXY*{ IPRINT-1) 

DO 25 TX=^1,NX 
WRITE (6, 300) 

00 25 IW=l,NW 
I WX=IW*IX 
INO=IND+l 

WRITF(6,200) (SOURCE (KK) ,KK=IND, INDFINtNWX) 

CONTINUE 
FQRMAT(8E15.5 ) 

FORMAT!/) 

RETURN 

CONTINUE 

RETURN 

CONTINUE 

RETURN 

CONTINUE 

WRITE(6,770) 

FORMAT(///'D'lSTRIBUTION OF AA(I,J)’/) 

DO 77 1=1, NTOTAL 
WRITE(6,771) 1 
N1=I 

^42 = NTOTAL*NTnTAL 
FORMA T(2X, ’ IN'OEX=* , 12 ) 

IF(NY,LE,4.0R,NY,GE.B ) write (6, 772 ) ( AA(K) , K=N 1 , N2 , NTOTAL) 

FORMAT! 8E1 5.6/3 El 5,6/ 8F 15 «6/ 8 FV5, 6/aE 15.6 ) 

IF(NY.E0.5)WRITE{'6',77 5) ( AA(K) ,K=NI, N2, NTOTAL) 

F0RMAT(5E15.6) 

IF{NY.G0.6) WRITE(6,776) ( AA ( K ), K=NI , N2 , NTOTAL > 

FGRMAT(6E15.6) 

IF (NY. EO. 7) WRITE! 6,777) (AA(K) ,K=N1,N2, NTOTAL) 

F0R’^AT(7E15. 6) 

CONTINUE reproducibility of the 

ORIGINAL PAGE IS POOR 
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return 

5 CONTINUE 
RETURN 

6 CONTINUE 
,WRITE{6t88l) 

881 FORMAT! ///2X» »THE DISTRIBUTION OF SURFACE NORMAL*/) 
NXW=NX*NW 
NXWY = NXW=!=NY 
DO 883 IX=1,NXW 

WRITE (6, 882 ) ( SOURCE tKK ) , KK= I X ,NXWY ,.NX-K) 

832 FORMAT(8E15,6) 

883 CONTINUE 
RETURN 

7 CONTINUE 
RETURN 
END 
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subroutine CFBUG(K) 

WRITE {6t DK 

I F0RMAT(2Xt 'ERROR C0nE=*,I2) 
RETURN 
END 
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C 


C 

c 


2 


c 


SUBRnUTIME VPLKf' 

C OMMON/ ?: Z Z,l / iMX , N!Y , N Z , N W » R EF L E N t SP AN , K SY M MY , KSYMM Z , NS Y MMY , NS YMMZk 
C 0 MM 0 M/ZZZ 2 /TAUt ALFAt TANGLE i T ANGTE , CHORD, NTOT aL t UMACH 
CCMM 0 N/ZZZ 6 /XPC (2 50 '), YPC( 250 ) ,ZPC( 250 ) 

COMMON/ZZZ 7 /XP 1 ( 250 ),YD 1 ( 250 ) ,ZPU 250 ) , XP 2 ( 250 ) ,YP 2 { 250 ) , 

IZP2(2 50) ,XP3(250) ,yP3>(-250)aP3(.250) 

C 0 MM, 0 N/ZZZ 8 /AA{ 2500 ) , SOURCE ( 250 ), S IN ABC , COS ABC, -AL F ABC 
.C 0 MMnN/ZZZ 9 /XPP (2 50 ), YPP( 250 ) , ZPP (2 50 ) , XPM( 250 ) ,YPM‘( 250 ) 

I, ZPM( 25 0 ) ,XNP (250 } , YMP ( 250 ) , ZMP ( 2 50 ) , XHM( 250 ) ,YMM( 250 ) , 

1ZMM{250) , I WAKE (250) 

COMMON/ ZZZ 10 /JNXB ( 250 ) tNXWAKE ,WAKEIN 

COMMON/ZZZll/VHKX (2500) , VHKY ( 2 500 ) , VHK Z'{ 2500 ) , VKX'(.250 ) ,VKYX250) , 
1VKZ( 250) 

C 0 ,MM 0 N/ZZZ 12 /VXWAKF (250 ) , VYWAKE (2 501 , VZV!AKE(-. 250 .). 

COMMON/COKTR/ NITER 

DOTPRG(Xl, Yl,Zi ,X 2 ,Y 2 ,Z 2 )=Xl*X 2 4 -Yl*Y. 2 -t-Zl>i'Z 2 

PRQMIX(XX 1 ,YY 1 , ZZ 1 ,XX 2 ,YY 2 ,'ZZ 2 ,XX 3 , YY 3 , ZZ 3 ) = (-YY 2 ’!=ZZ 3 ~YY 3 X'ZZ 2 )*XX 1 
l-( XX 2 =^ZZ 3 -XX 3 YZZ 2 )*YY 1 -^(XX 2 =^'YY 3 -XX 3 =?=YY 2 )*ZZ 1 
NT 2 S=NTOTAL’ 5‘>!=2 

NYP=MY +1 
NWT^NXWAKE*NY 
NWTP = NWT 4 -NXWAKE 
NTBW=NTOTAL + NV.T 


DO 2 1=1,250 
VXWAKE( I )=0c 
VYWAKE( I )=0. 

VZWAKE( n=o. 

CQNTTNUE 

DO 250 JNXBW=1,NTBW 

DO 250 INXW=l,NWr 

DO 250 ISYMMY=1 ,MSYNMY 

,D0 250 rSYMMZ=l ,NSY^^MZ 

SI‘GNY=3--2*ISYMMY 

SIGNZ=3.-2*ISYMM7 

JNXW=JNXSW-NTCTAL 

IF('JNXBW.LE,NT0TAL) JMX=JNX8W 

I'E (JNXBW.GT. NTOTAL )JNX=JNXBTJNXW) 


INDEX=='NTGTAL + TNXW 
01X=XPP( JNXB-W)-XMM (INDEX) 
qIY^YPP ( JNX8W)-YMM ( index ) *S IGNY 
01Z=ZPP( JNXBW)-ZMM( tndEX-)*SIGNZ 
02X = XMP( JNX8W)-XYP(TnceX > 

Q2Y=YMP( JNXBU)-YMM( INDEX) «SIGNY 
Q2Z=ZMP( JNXew)-ZYp( index ) ^SIGNZ 
Q3X=XMM{ JNXBW) -XMW (INDEX) . 

Q3Y = YiMM( JNX3W)-YMM{ INDEX )*SIGNY 

•03Z'=zmm.('Jnxbw)-Z‘'^m( i'ndex)*s IGNZ 

QA-X=XPM( JNXBW )-XMM ( INDEX) 

D4Y=YPM ( JNXBW ) -YMw { INDEX ) *S IGNY 
Q4Z= ZPM ( JNX RW ) -Z X P (I N DEX )'*S IGNZ 
QlQl=DOTPRn(QlX ,Q I Y , Q IZ , C IX , G lY , Q 1 Z ) 
Q2Q2=OOTPRO(C2XiC2Y,Q2Z,C2X,Q2Y,Q2Z) 
Q3Q3=D0TPRG(Q3X,Q3Y,Q3Z, Q3XrQ3Y,03Z) 
Q4Q4=nQTPRG(G4X,Q4Y,Q4Z,Q4X,Q4Y,Q4Z) 



- 89 - 


C 


250 


A 


3 


753 

5 


QrQ2=DOTP.R0{Ql'X,mYtOlZ,C2X,O2Y,O2Z) 

QlQA=0aTPR0 (C1X»01Y,Q1Z ,G4X,C4Y,04Z) 
Q2Q3=DOTPRa(Q2X,O2Y,0'2ZTQ3X,Q3Y,Q3Z) 
Q30A=0aTPRQ{03X,Q3Y,03Z,QAX,QAY,04Z) 

01=$QRT(qiGl) 

Q2=S0RT(Q2Q2) 

03=S0RT(Q3Q3) 

Q4=SQRT{04G4) 

Q4iX^04Y*01Z-Q4Z«QlY 
041Y=-(Q4X=<=G17-C4Z#CIX) 
q41Z = 04X=^^01Y-C4Y*0lX 

Q41SQ=nOTPROlQ41X,q41Yt041Z,O'^lX,Q41Y,Q41Z) 

0 1 2 Q 1 Y =>= Q2 Z - C 1 Z'-i^C2 Y 
0 1 2Y=- { Q 1 X’i'O 2 Z-Q IZ’^'GZ X ) 

Ql2Z=qiX=f'Q2Y-QlY«02X 

Q12SQ=nOTPRQ(Ql2Xt012Y,012Z,Q12X,Q12Y,Ql2Z) 

Q23X=Q2Y*03Z-C2Z*03Y 

023Y=-(02X’i'03Z-C2Z’5‘G3X) 

Q23Z=Q2X=:'q3Y-G2Y*Q3X 

Q23SQ=nnTPR0C C23Xt023Y,Q2 3Z,Q23X,Q23YTQ23Z) 

Q34X=03Y=!=C4Z-C3Z*C4Y 

Q34Y=-(Q3X*04Z-C3Z^Q4X) 

034Z=03X=i'Q4Y-C3YY04X 

Q34S0=DGTPRG(Q34X,Q34YiQ34Z,Q34X, Q34Y,Q34Z) 

PARTI=0* 

IF'(Q41SQ.NE.O. ) PARTl= U Q4Q4-Q1Q4 ) /q44- { Q l.Ql-0104 3 / Q 1 ) /041SO . 
P ART2=0. 

IF(012SO.NE.O,) PART2= ( (.QlCl-Ql Q2 ) /Q 1+ ( Q 2Q-2-Q IQ2 ) / Q 2 3 / 0 12 SO 
PART3=0. 

IF(Q23S0-NE, 0, ) PART3= {( 02C2-Q2Q3 ) /02 «■ { Q303-Q2Q3 } /Q3 ) /023 SQ 
PART4=0* 

IF{Q34SO.NE.O. 3 PART4={ ( 03 03-0304) /Q3+ ( Q4Q4-C3C4 ) /C4 )/ q34SC 

VX=Q4iX*PARTl ^012X*PART2^■Q23X^PART3f0 34X*PART4 

VY= CQ41Y=XPARTU-qi2Y’l'PART2 + 023Y=f‘PART3+Q34Y’?'PART4}«SIG''lY 

VZ = q41Z*PARTU-Q 12Z^PART2+Q23Z^PiRT3i-g34Z*PART4 

VXWAKFt INXW3^VXWA!<E (INXW)+VX=!‘SCURCF(JNX3 

VYWAKE( I\'XW)=VYWAKE( INXW)+VY^SOURCE( JNX) 

VZWAKEC INXW)=VZWAKE { I NXW 3 +VZ-S0URC6 ( JNX 3 
CONTINUF 
CALL VELPP 
DG 4 I=1,NTGTAL 
VKX{ I 3=0. • 

VKY( I )-=0. 

VKZ( I 3=0. 

CONTINUE 
DO 3 T=l,NTOTAL 
DO 3 J=1,NT0TAL 
NNN = l4-{ J-ll'^NTCTAL 
VKXn 3 = VKX{ 1 3 + SOURCE (J)*VHKX{NNN) 

VKYt I)=VKYn ) <-SCURCF{ J)*VHKY(NNN3 
VKZn 3=VKZ( I ) +SCURCE( J 3*VHKZ { NNN) 

CONTINUE 

IF (NITER. EQ. 13GG TO 753 
IF(NITFR.LE.103G0 TO 2000 
WR IT E { 6 T 5 3 

>"ORMAT(/10Xt ' this IS THF X-WING VELOCITY’/) 

CALL PRINTVI Vi<XfNX,NY) 


rEPRODUCIBIL^ OF 
ORIGINAL PAGE IS POOR 
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WRITF(6t6) 

6 FORMAT(/IOX, ’THIS IS THE X-WAKE 
CALL PRINTVIVXWAKE, NXWAKE,NYP) 
WRITE (6t 7) 

7 F0RMAT(/10X, • THIS IS THEV-WING 
CALL PRINTVI VKY,NXt NY) . 
WRITE(6,8) 

8 ■ FORMAT {/lOX, ’ THIS IS THE Y-WAKE 

CALL PRINTVI VYWAKE,NXWAKE,NYP) 
WRITE (6, 9)' 

9 T^ORMATI/lOX, 'THIS IS THE Z-W I HG 
CALL PRINTVI VKZ,NfX,NY) 

WRITEI6, 10) 

10 FORMAT(/10X, * THIS IS THE Z-WAKE 
CALL PR I N TV ( V ZW A KE t NX WA KE , NY P ) 

2000 CONTINUE 
RETURN 
END 


VELOCITY*/) 

VELOCITY*/) 

VELOCITY'*/) 

•V*EL0C1T-Y'-/) 

VELCCTTY'/) 
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SUSROUTI’NE PR INTV (VECTOR, N1 , N2 ) 

DIMENSION VECTOR! 1) 

WRITE(6,3 } 

DO I IX=l,Nl 

WRITE (6,2) (VECTOR (I X+Nl* (I Y-1.) ) , IY=L,N2> 

1 CONTINUE 

2 F0RHAT(8E15,6) 

WR.ITE(6,3) 

3 FORMAT!/) 

RETURN 

END 
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S’UBRnUTINE COEFF 

C0MMQN/ZZZ1/MX,NY ,NZ, NVJ,REFLEN,SPAN^KSYMMY»KSYMM2tNSYMMY,NSYMM^ 
C0MM0N/ZZZ2/TAU» ALFA, TANGLE, TANGTF, CHORD, NTOT AL , UMACH 
CQMMON/ZZZ6/XPC(25O),YPCf250) ,ZPC1250) 

C0MM0N/ZZZ7/XP1 C 250), YP 1(250) , ZPl ( 250 ) , X°2 ( 250) ,YP2(250) , 
l.ZP2(250) ,XP3( 250) ,YP3‘(250) ,ZP3(250) 

C0MMQN/ZZZ8/AA{2500) ,S0URCE(250) , S I N ABC , COS ABC, ALFABC 
COMHON/ZZZ<5/XPP‘{2 50), YPP (250) , ZPP ( 250 ) , XPM ( 250 > ,YPM(250) 

1 , ZP M ( 2 SO ) , XMO ( 2 50 } , YM P ( 2 5 C ) , ZMP ( 2 5 0 ) , X MM ( 2 50 ) , YMM ( 25 0 ) , 

IZMM{250) , Iv'JAKE{250) 

COMMON/ZZZ10/JNXB(250 ),NXWAKF,WAKE1N 

COM-^ON/ZZZ11/VHKX(2 500) ,VHKY(2500 ), VHKZ(2500 ) ,VKX { 250 ) , VKY ( 250) , 
-LV'I^Z (250) 

OT'HENSION XUN0RM(250) , YUNGRM { 250) ZUNORM { 250) 
■D‘0TPR0(X1,YL,Z1,X2,Y2,Z2)=XI*X2+Y1=>‘Y2J:Z14=Z2 

PR0MIX(XX1,YY1, ZZl, XX2, YY2, ZZ2, XX3, YY3 , ZZ3) = ( YY2*Z23-YY3*ZZ2 1 «XX;. 
I--(XX2*ZZ3-XX3*ZZ2 )*YY14-(XX2=?^YY3-XX3*YY2)*ZZ1 
NT2S=NT0TAL=?=>!^2 


00. 9 NNN=l,NT2S 
VHKX(NNN)=0. 

VHKY(NNN) -=0. 

VHKZ(NMN)=0. 

AA(MNN)=0« 

CALCULATION OF THE SURFACE NORMAL 

i 

DO 140 JNX-1,NTGTAL 
C 

XD1=XPP{ JNX)-XMM( JNX) 

YD1=YPP( JNX)-YMM(JMX) 

ZD1=ZPP( JNX)-ZMY(jmx) 

XD2==XMP( JNX)-XPM{ JNX) 

YD2=YMP( JNX)-YPy{ JNX) 

ZD2=ZMP( JNX)-ZPW( JNX) 

CRX=Y01*ZC2-ZDl=frYD2 
CRY=- {X0l*ZC2-ZDl^Xn2 ) 

CRZ = XDl=i'YD2-YDl*XD2 

ABN=S0RT (DOTPRO (C RX , CRY , CRZ , C RX ,CRY , CRZ ) ) 
XUNORM{JNX)=CRX/ABN 
YUNQRM(JNX)=CRY/A6N 
ZUNORM(JNX)=CRZ/A8N 
1-40 CONTINUE 

NTBW^NTOTALfNXWAKe^NY 
DO 250 JNX8W = 1, N'TBW 
C 

DO 250 INX=1,NTCTAL 
DO 250 ISYMMY=1,NSYMMY 
00 250 ISYMW,2=i,NSYMMZ 
STGNY=3.-2=5'ISYMVY 
SIGNZ=3*-2’!'ISYMw2 
C 

JNXW=JNXBW-NTOTAL 
IF( JNXfiW.LE.MTDTAL) JNX=JNXBW 
I F( JNXBW.GT.NTDTAL) JNX= JNXB ( JNXW) 
NNN=INXi-(JNX-l ) *N TOTAL 
C 

Q1X=XPP{ JNXBW)-XPC( TNX) 
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0-lY=YPP'(>JNXBW)~YPC( INX)*SIGNY 
Q1'Z=ZPP C JNXBW)-ZPC MMX)’!'SIGMZ 
Q2X=XMP( JNX?vn-Xt>C( INX) 

Q'2Y=YMP{ JNXBW)-YPCnNX)»=SIGNY 
02Z=ZMP( JNX8W )-ZPC (INX)^SIGNZ 
Q3X=XMM( JNX6W)-XPCUNX) 

Q3Y=YMM{ JNXBW)-YPC( INX)=^=SIGNY 
Q3Z = ZMM( JNXeW)~ZPCnNX) *SIGNZ 
Q4X=XPM{ JNXBW)-XPC( IMX) 

Q4Y^YPM{ JNX8W )-YPC ( TNX^ =^=SIGNY 

Q4 Z= ZPM ( J NX 8W ) - Z P C (IN X ) « S I GN Z 

Q101=OOTPRO(Q1X,Q1Y,01Z»Q1X,G1Y,Q1Z) 

O202=DOTPR0(Q2X,Q2Y,Q2Z,Q2XtQ2Yt02Z) 

Q3Q3=n0TPR0(Q3X,Q3Y,Q3Z,Q3X,Q3Y', 03ZJ 

0404=n0TPR0 ( Q4X ,Q4Y,Q4Z ,.04 Xt C4Y404Z ) 

0102=nOTPRC(01X,QlYTQlZTG2X,G2Y,Q2Z) 

QlQ4=DOTPRG(Q.iXTQlY,OlZtG4XtQ4YyQ4Z) 

0203=D0TPR0(Q2XTQ2Yt02ZtQ3XtQ3Y,Q3Z) 

Q3Q4=OOTPRQ(G3X,G3Y,Q3Z,G4X,Q4Y,Q4Z) 

Ol^SQRT(OlQl) 

02=SQRT(02G2) 

03=S0RT-(Q3Q3) 

Q4=S0RT(Q4Q4) 

Q41X=Q4Y*QlZ-C4Z’t‘QlY 

Q41Y=-(04X*G1Z-C4Z*Q1X) 

Q4iZ = 04X=<=01Y-G4Y*QlX 

Q4lSO=nOTPRO(C41X ,04LYt Q41Zf Q41X, Q41Y,Q^1Z) 

Q12X=QL-Y*02Z-Q1Z*G2Y 

Ql2Y=-(OlX=i=02Z-GlZ*Q2X) 

Q12Z=Q1X^'Q2Y-G1Y*Q2X 

Q12SO=DOTPRO( QL2X ,012Y, Q12Z,Q12X, C12YtG12Z) 

Q23X=Q2Y*03Z-C2 2^«Q3Y 
023Y=-(02X^Q3Z-C2Z*Q3XI 
023Z = Q2X>!=Q3Y-g2Y^<03X 

023SQ = COTPROCG23X,Q23Y, Q23Z,Q23X-, Q23YtQ23Z) 

Q34X=Q3Y*C4Z-C3Z*Q4Y 

Q34Y=-(03X*Q4Z-C3Z>f=04X) 

034Z-Q3X’;'G4Y-G3Y=«=G4X 

Q34SQ=DOTPRO(Q34X ,Q34Y , Q34Z , Q34X , Q34Y , Q34Z) 

C 

PART1=0. 

IR{041SQ.NE.O. )PART1= { (Q4G4-QIQ4) /Q4+ { Q1 Q 1-01Q4 ) /Q 1 ) /-Q41 SQ 
PART2=0. 

1F(Q12SQ.NE.0. ) PART2= { ( Q1 C1-G1Q2 ) /0l + ( Q2Q2-Q1Q2 ) / Q2 ) /012S0 
P ART3^0. 

IF(02 3SQ.NE.O, ) PART3= ( (-02C2-G2Q3 ) /Q2-^■{Q3Q3-Q2Q3 ) /Q3 ) /Q23SQ 
PART4=0, 

IF(Q34SQ.NE.O.) PART4= ( ( 03Q3-0 3Q4 ) /O 3 K Q4Q4-Q3Q4 ) /C4 ) / 034SG 
VX=Q4lX*PARTl+Q.12X-PART2■^C23X*oART3^■Q34X*PART4 
VY=(041Y*PARTUQl2YS‘PART2 + Q23Y*PART3i-Q34Y*PART4)*SIGNY 
VZ^Q4LZ*PARTl^■Ql2Z=^'PART2^G23Z=^PART3^034Z*PART4 
C VHKX(NNM)=VHKX(NNN) +VX 

C VHKY(NNN)=VHKY{NNN)+VY 

VHKZ{NNN)=VHKZ(NNN)fVZ 

FACTnR = DOTPRG(VX,VY,VZ,XUNGRM( INX) ^ YUNORM ( I NX ) , ZUNORN ( INX) ) 
AA'{NNN)^AA (NNN } +F ACTOR 

THE NEXT FEW LINES ARE A FEW CHECK STATEMENTS 
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C SOURCE ( INX)=SOURCE-{ INXX-FACTOR 

250 ■ CONTINUE 
C . WRITE (6, 152) 

152 F0RMAT(20X, ' THIS IS THE MATRIX AAV/) 

C WRITE(6, 15U(AA(I),I=1,NT2S) 

C WRITE(6,153) 

153 FnRMATt//20X,,’THIS IS THE MATRIX SOURCE’’//) 

C WRITF(6,151) (SOURCE (I 5 , I^ 1 , NTOT AL ) 

151 FORMAT(8F15.6) 

no 154 T=1,NT0TAL 

SOURCE { I )=-{XUNORM ( I ) =«'CGS-A!3C-s- ZUNORM ( I ) ^SINABC ) 

154 CONTINUE 
•RETURN 
ENn 
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SURROUTIMF V*=LPP 

COMMON/ ZZZI/NXt NY »MZ, NW , RpPLEN, SP AN , KSYMMY» KSYMMZ »NS YMMYt NSYM m^ 
COMMON/ZZZZ/T^U, ALFA, TANGLE, TANGTE, CHORD, NTOTALtUM ACM 
C0MM0N/ZZZ6/XPC(2 50 >, YPC( 250) ,ZPC (250) 

COMMON/ ZZZ7/XP1 (250), YPl (250) , ZPl i 250 ) , XP2 ( 250 ) , YP2 ( 250 ) , 
1ZP2(250) f XP3( 250) ,YP3 (250)> ZP3(250) 

CDHMON/ZZZa/AA( 2500), SOURCE (250), SIN ABC , COS ABC, ALF ABC 
C0MM0N/ZZZ9/XPP (250), YPP (250) ,ZPP (250 ) , XPM ( 250 ) , Y PM { 259 ) 
1,ZPM(250) ,X.MP(250},Y-MP( 250) ,ZMP(250) ,XWM( 250) ,YMM(250) , 

1ZM.M(250) , I WAKE (250) 

CGMM0N/ZZZ10/JNX3(250) ,NXWAKE,WAKEIN 

C0MM0N/ZZZll/VHKX(250G) ,VHKY( 2500),VHKZ(2500) ,VKX(250) ,VKY(250) , 
1VKZ(250) 

C0MMON/ZZZ12/VXWAKEC250) , VYWAKE(250) ,VZ WAKE (250) 
DOTPRO.{X1,Y1,-71,X2,Y2,Z2 )=X 1*X2 4-Y 1 *Y2 fZ I*Z2 

PROMTXI-XXl ,YY1,7Z1, XX2, YY2,ZZ2,XX3, YY3, ZZ3) = ( YY2=f=.ZZ3-YY3’!'ZZ2)«xxi 
l-(XX2*ZZ3-XX3’!<ZZ2)=<=YYU-( XX2^YY3-XX3’i=YY2)>l'ZZl 
N‘T2S=NTOT'AL=4‘*2 

N.WT=NX.WAKE*NY 

NVITP1=NWT^-1 

nwtp^nwt+nxwake 

NTBW=NTOTAL«-NWT 

INITIALIZATION OF THE WAKE VELOCITY AT THE CORNERS 

DO 250 JNX8W=1,NTEW 

DO 250 INXW=NWTP1 ,NWTP 

no 250 ISYMMY=l,NSYNMY 

DO 250 ISYMMZ=1 ,NSY«MZ 

SIGNY=3.-2*TSYMKY 

SIGNZ=3.-2«ISYMMZ 

JNXW= JNXSW-NTGTAL 

IF( JNXBW.LE.NTOTAL) JNX=JNXBW 

rF{ JNXBW.GT.NTDTAL) JNX=JNXB(JNXW) 

INOEX=NTOTALi-INXW~NXWAKE 
Q1X=XPP< JNX8W]-XMP( INDEX) 

0LY=YPP-( JMXBW]-Y^P( INDEX)=?SIGNY 
QIZ^ZPP ( JNXBW ) -ZMP ( INDEX ) IGNZ 
02X=XMP( JNXPW)-XMP( INDEX) 

02Y=YMP{ JNXBW )-YMP( INDEX )*SIGNY 
02 Z=ZMP (JNXBW) “ZMP ( INDEX ) >^S IGNZ 
Q3X=XMM( JNXBW )-X^'P( INDEX) 

Q3Y=YMM{ JNXBW )-YMP (INDEX ) *SIGNY 
Q3Z=ZMM { JNXBW) -ZMP ( INDEX ) =*S IGNZ 
Q4X=ypM{ JNXBW) -XMP( INDEX) 

Q4Y=yPM( JNXBW )-YMP( IND E X ) »^ S I GN Y 
04Z=ZPM{ JNXBW)-ZMP( INDEX) =<'S IGNZ 
Q101 = DQTPRO(Q1X,Q1Y,01Z,C1X,G1.Y,01Z) 

0202= DOT PRO ( Q2X, 02Y, 02Z, Q2X,Q2Y,Q2Z) 

0303=D0TPR0(Q3X,Q3Y,Q3Z, C3X,Q3Y,03Z ) 
Q4Q4=D0TPR0(Q4X,04Y,Q4Z,CAX,G^Y,Q4Z) 
Q102=D0TPRa{QlX,giY,QlZ,G2X,Q2Y,Q2Z) 
Q104=DnTPRQ(giX,01Y,01Z,G4X,Q4Y,04Z ) 
02Q3=0OTPRn(02X,Q2Y,Q2Z,03X,Q3Y,.O3Z) 

03Q4=nOTPRn(C3X,03Y,Q3Z, G4X, Q4Y,Q4Z) 

01=S0RT(Q1Q1) 
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Q2=S0RT{Q2Q2) 

03=SQRT (Q3Q3) 

04=SQRT(0-^Q4) 

041X=Q4Y=?QlZ-Q4Z=i'01Y 

041Y=-(Q4X*QIZ-Q4Z*Q1X) 

04iZ=Q4X*QlY-C4Y«QiX 

041SO=noTPR0(Q41X,Q4lY,041Z,Q41X, Q41Y,Q41Z) 

012X-Q1Y«Q2Z-C1Z*C2Y 

Q12y=-(01x-Q?Z-QLZ*Q2X) 

012Z=01X«Q2Y-C1Y*Q2X 

012SQ=nOTPRO{G12X»C:i2Y, Q12 ZtQL 2X, Q12YtQl'2-Z ) 

ORSX^QZY^QIZ-QZZ^J'GSY 

Q23Y=-(02X'-l'Q3Z-C2Z*Q3X) 

023Z=02XYQ3Y-G2Y=<=03X ‘ 

Q23S0=nOTPRO(Q23X,Q23Y,g23ZTQ23X,Q23Y,C23Z) 

Q34X=Q3Y*04Z-G3Z’-C4Y 

Q34Y=-(03X=i=Q4Z-Q3Z=!' G4X) 

034Z=03X*04Y-Q3Y*Q4X 

Q34SQ=nOTPRO(C34X,G34YT034Z,Q34X, G34YyQ34Z3 
C 

parting. 

TF(Q41SO.NF.O. ) PABT1= ( (Q4Q4-Q1Q4) / Q4+ ( 0 10 1-Q 1Q4 ) /Q 1 ) /041 SO 
PART2=0.‘ 

IF ( 012S0.NE.O. ) PART2= I { Q1 C 1-Q 1Q2 ) /Q 1 K Q 2Q2^Q IQ2 ) / 021/ 0 12 SO 
PART3 = 0. 

IF(Q23SQ.NE.0. )PART3= ( ( Q2G2-Q2Q3 ) /02i- { G3Q3-Q2G3 ) 7 Q3 ) /Q23SQ 
PART4=0. 

IF (034S0.NF. 0. ) PART4= ( ( G3C3-Q3Q43 /Q3 + I 04Q4-Q3Q4 ) / Q4 )/Q34SQ 
VX=Q41X^PARTl+0 12X*PART2^-Q23X*PART3+Q34X*PART4 
VY={Q41Y*PART1 + G12Y*P ART2+Q23Y#PART3+-Q34Y=!'PART4) *S IGNY 
VZ = Q41Z*PART1 + 012Z*PART2+-C23Z*PART3 + Q34Z>!'PART4 
IF( INXW.LE.NWDGO TG 250 
VXWAKEI INXWI^VXWAKEI INXW) fVX’^'SCURCEI JNX) 

VYWAKEI INX'W) = VYWAKEUNXW) +VY*SCURCE( JNX) 

VZWAKEI INXW)=VZWAKE(INXW) fVZ*SOURCE(JNX) 

250 CONTINUE 
RETURN 
END 
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subroutine ITER 

COMMON/Z7Zl/NX,NY,NZ,NWtR'"FLFN,SPAN,KSYMMY,KSYMMZ,NSYMMY’, NSYMMZ 
COMMON/ Z Z Z2 / T AU t ALFA, TANGLE, TANGTEt CHOP D , NTDTAL , UMACH 
COMMON/ZZ26/XPC(250), YPC(250) ,ZPC(250) 

C0MM0N/ZZZ7/XPU250), YPK250) , ZPH 250 ) , XP2 ( 250 ) ,YP2'(250) , 
1ZP2I250) ,XP3(250I ,YP3(250) ,ZP3(250) 

COMMnN/ZZ28/AA(900) ,SOURCE(250) , S I N ABC , COS A8 C , AL F A BC 
COMMON/ ZZZ9/XPP (250) , YPP( 250) , ZPP ( 2 50 ) , X°M( 2 50) ,YPM(2'50) 

I, ZPM(250) ,XMP (250) , YMP(25C) , ZMP (250 ) , X«M { 250 ) , YMM { 250 ) , 
1ZMM(250) ,IWAKE(250) 

CGMM0N/7.ZZ10/ JNX8(250 ) , NXW AK E , W AK E IN 

C0MM0N/ZZZ11/VHKX(900) , VHKY ( 900 ) , VHKZ { 900 ) , VKX ( 250 ) ,VKY( 250) , 
1VKZ( 250) 

COMMON/ ZZZ12/VXWAKE (250) ,VYWAK,F (.2 5.0.)., yXW.AK.E,(.250.0. 
COMMON/CQNTR/NITER 

DIMENSION XXX(250) ,YYY(250) ,ZZZ(250) ,VXW(250) , VYW ( 250 ) , VZW ( 250 < 
DIMENSION INDICA(IOO) 

nOTPRO(Xl,Yl ,Zl,X2,Y2,Z2)=Xl*X2+Yl’f^Y2FZl*Z2 


NXWAKP==NXWAKE + 1 

ALFAR=ALFA*3.14159/180. 

TANALF=TAN(ALFAR) 

NYP=NYH 

00 1 IX=1,NXWA'KE 
DO 1 IY=1,NY 
IXP=IXi-l 
TYP=IY-+1 

IELEM=NT0TAL+IX4-( IY-1 )*NXWAKE 
lNODEl=IXP+( !YP-1)*NXWAKP- 
IN0nE2=IX +{ IYP-1 )*NXWAKP 
IN0DE3=TX <-(IY -1)*NXWAKP 
INQDEA^IXPK lY -D’i'NXWAKP 
C 

XXX( INODEU^XPPdELEM) 

YYY( INODEL)=YPP (I ELEM) 
2ZZ(IN0DE1)=ZPP ( IELEM) 

C 

XXX( IN00E2) = XMP (IELEM ) 

YYY( IN0DE2)=YMP(IELEM) 

ZZZ( IN00E2)=ZVP( IELEM) 

C 

XXX( IM0nE3)=Xwy.nELEM') 

YYY( IN0nE3)=YMM( IELEM) 
ZZZ(IN0DE3) = ZMK( IELEM ) 

C 

XXX ( IN00E4)=XPM( IELEM ) 

YYY { I N0DE4) =YPM ( IELEM ) 

ZZZ( IN0DE4)=ZPM(IELEM) 

C 

1 CONTINUE 

DO 50 IX=1,NXWAKP 
00 50 1Y=1,NYP 
INniCA(TY)=IY*NXWAKP 
IN0DE=IX4-(IY-1)*MXWAKP 
INDEXl=INQD6-( lY-1) 

VXW( I NODE) =VXWAKE( INDEX 1 ) 

VYW( INODE)=VYWAKE( INDEXl) 
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50 
C 
C 

c 

51 
C 

52 


5.3 


100 

201 

202 

203 

204 

205 

1000 

C 

C 


C 

C 

C 

2000 


VZwnN006)=V7WAKFnNnEXl) 

IF( IX.EO.NXWAKD ) V XW { I NODE ) = 0 • 
IP( IX, EQ.NXWAKP)VYW{ INn0E}=0. 
IF( IX.EQ.NXKAKP ) V ZW { T NODE ) =0 . 
IF( IX.EQ.l )VZW(!NjnOE)=-TAKALF 
IFCIX.EQ, 1} VYW(INGOE) =0. 

I F{ IX ,E0. 1 )VXW( I^1GDE) = 1, 
CONTINUE 


WRITE (6,51) 

F0RNAT(/2X, ’ PRINTOUT OF THE 
CALL PRINTV( VXW,MXWAKPtNYP) 
WRI-TF(6,52) 

F0RMAT(/3X, ' PRINTOUT OF THE 
CALL PRINTV(VYW,MXWAKP,NYP) 
■WRITE (6, 53) 

F.ORNATI /3X, ' PRI-NTOUT OF THF 
CALL PRINTV{VZW,NXWAKP,NYP} 
IF ( NITER, C-T. DGG TO 1000 
WRITFI6, 100) 

FORMAT (/3Xt ’ PRINTOUT OF THE 
WRITE(6,20i) 

FCRMAT{3Xt ’ ITERATION IN THE 
CALL PRINTV(XXX,NXWAKP,NYP) 
WRITEI6, 202) 

FORMAT (/ 3 Xt ’ PRINTOUT OF THE 
WRITF(6,203) 

F0RMATC3X, ‘ITERATION IN THE 
CALL PRINTV(YYYtMXWAKPtNYP) 
WRITE [6t204) 

FORMAT! /3X, ' PRINTOUT OF THE 
WRITE (6t205) 

F0RMAT(3Xt • ITERATION IN THE 
CALL PRINTV(ZZZtNXWAKP,NYP) 
CONTINUE 


WAKE X-VELOCTTYV) 

WAKE Y-VELOCITY‘/ ) 

WAKE Z-VELOCTTY'V) 

WAKE CORNER COORDINATES BEFORE') 
X-DIRECTTON*/} 

WAKE CORNER COORDINATES BEFORE') 
Y-DTRECTICN'/ ) 

WAKE CORNER COORDINATES BEFORE’) 
Z-DIRECTION •/) 


DO 3 IX=1,NXWAKP 
DO 3 IY=1,NYP 
INODE=IX+( IY-1 )*NXWAKP 
IFnX.EO.NXWAKP)GG TO 3 
R=WAKEIN 

VELTnT = SQRT(noTPRC{ VXWl'lNODE) ,VYW{INODE) tVZWC INODEJt 
IVXW(INO nF) ,VYVH INODE) ,VZW (INODE) ) ) 

IF( VELTOT.EQ.Q. )CALL CEBUG(50) 

DELX^RX:VXW{ INODE) /VELTOT 
DELY=R*VYW ( I NODE) /VFLTQT 
DELZ=R=!<VZW( INQOE) / VELTOT 

INQDPl=INOnE+l 

I F( INOOPl.EQ. INDtCAI lY) )GC TO 2001 

XXXnNODPl} = XXXtINODPl) ^OELX 

YYY(INOnPl)-=YYY(INOnE ) +DELY 

TFIIX .EO.l.AND, IY.EQ.NYP)GO TO 2000 

Z7Z{ INODPD^ZZZUNODE) fOELZ 

GO TO 3 

CONTINUE 
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IF{NITER,P0.1 )ZZZnNODPl)=ZZZnNOnE) fDE.L2 

c • • 

C BRING THE WHOLE VORTEX IN LINE WITH THE LAST Z 

C 

INOEXl-INODPl+1 
INDFX2==NXWAKE-M lY-1 I^i'NXWAKE 
C on 2 IMOEX^IN0EXl,TNnEX2 

C2 ZZZ( INnEX)=2ZZnNODPl) 

C 

GO TO 3 

2001- CONTINUE 

YYY( JN0DP1)=YYY n NODE) 

■ ZZZUNOOoi ) = ZZZ (INODE ) 

C 

3 CONTINUE 
C 

DO 4 IX=1 ,NXWAKE . 

DO 4 IY=1,NY 

IELEM = NTOTAL+IXt( IY-1 )=^NXWAKE 

IXP=IX+1 

IYP=IY+1 

IN0DEl = IXPJ-( lYp-l )YNfXWAKP 
IN0DE2=IX i-( I YP-1 )*NXWAKP 
IN00E3=IX +(!Y -D^NXWAKP 
IN0DE4=TXPM lY -1)=«'NXWAKP 
C 

XPP{ IELEM3 = XXX( INODEl ) 

YPPd ELEM)=YYY{ T.NQ0ET )- 
ZPP(IELEM)=ZZZ( INODEl ) 

C 

XPM(IELEM)=XXX( INC0E4 )■ 

YPM(IELEM)=YYY{ INC0E4) 

ZPM(IELEM)=ZZZ{ IN0DE4) 

C 

IF{ IX.EO.l )G0 TO 4 
C 

XMM(I ELEM)=XXX( ING0E3 ) 

YMM{ IELEN)=YYY( TN0DE3 ) 

ZMM( I ELEM )~ZZZ ( INCDE3 ) 

C 

XMP{ TELEM)=XXX{ TN0DE2 ) 

YMPd ELEM)=YYY{ INODE? ) 

ZMP( IELEM)=ZZZ( INOOE2) 

C 

C 

4 CONTINUE 

IF (NITER. LF. lOlGO TO 738 
C WRITE{6,400)NITER 

400 F0RNAT{/3X, 'AFTERS l3>2X,dTERATI0NS» THEX-CORNER ») 

C WRITE{6r401) 

401 FORMAT (3X, ' COORDINATES OF THE WAKE ARE*/) 

C CALL PRINTV( XXX ,NXWAKP,NYP) 

WRITF(6,402)NITER 

402 F0RMAT(/3X, 'AFTER S I3t2X, « ITERATIONS, THE Y-CQRNER ' ) 
WRITF(6,403 ) 

403 FORMATOX, « CnOPHINATFS OF THE WAKE ARE'/) 

DO 601 TX=1,NXWAKP 

DO 601 IY=1,NYP 
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.lN0Df=lX4-( lY-D ’S'NXWAKP 

YYYnNnOF) = (l./(SPAN72. ) ) ^>YYY aNOCE ) 

601 CONTINUE 

CALL PRINTV{YYY,NXWAKP,NYP) 

WRITE{6,40A) NITE=^ 

404 ^0RMAT(/3X, ’AFTER ‘ tI3,2X, 'ITERATIONS, THE Z-CGRNER') 
WRITE(6t405) 

405 FORMATOX, ’COGRDIKATES OF THE WAKE ARE’/) 

CALL PR INTVf ZZZ ,NXWAKP,NYP) 

738 CONTINUF 
RETURN- 
END 
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■^3^ Pr$i%tou-fe of -Computer- Program ILSAT-7R 


SPECIFICATIGNS OF THE PROBLEM 

KiX= 7 
-NY^ 7 
NZ= 1 
NW= 1 

NTOTAL= 49 

KSYMMY= 1 
KSYHMZ= 0 

REFERENCE LENGTH^ l.OO 
SPAN/REF LENGTH = 8.00 

THICKNESS^ 0.0 

ALPA= 5.000 
ALFA3C= 0.0 

MACH NUMBER = 0.0 


NXWAKE= 11 ■ 

WAKEIN= 0.500 

_JTiiNGLEj= ^0..,0 

TANGTE= 0.0 


CH0RD= 1,00 



REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


PRINTOUT OF THE MAKE CORNER COORDINATES BEFORE 
ITERATION IN THF X-OIRECTION 


-0. 118755F-C6 
0.500000E 00 
0 . nonooE oi 

0.150000E 01 
0.200000E 01 
0.P500C0F 01 
O.BOOOnOE 01 
0.3E000OE 01 
O.AOOCOOF 01 
O.ASOOOOF 01 
0.500000F 01 
O.PBOOOOE 03 


-0.U8756F-06 
0.500000F 00 
O.lOOnOOF 01 
O-ISOOOOC 01 
o.2noonoE oi 
0.2’50000F 01 
0.3O0000E 01 
0. 350000= 01 
O.AOOOOOE 01 
0.A50000E 01 
0.500000E 01 
0.55000OF 03 


-0. 1187566-06 
0.500000E 00 
O.IOCOOOR 01 
0. I50000E 01 
0.2000006 01 
0.250000E 01 
0.300000F 01 
0.350000F Cl 
0.6000006 01 
0.650000F 01 
0.500000F 01 
0.5500C0E 03 


-0.118756F-06 
0.5000006 00 
0.100000= 01 
0.150000F 01 
0.2000006 01 
0.250000E 01 
0.3000006 01 
0.350000F 01 
0.4000006 01 
0.450000F 01 
0.500000= 01 
0.550000E 03 


0. U8756E- 

-06 

0.5000006 

00 

0. 100000 = 

01 

0. 150000 = 

01 

0.200000E 

01 

0.2500006 

01 

0. 3000006 

01 

0.350000= 

01 

0.409009E 

01 

0.450000F 

01 

0. 500000F 

01 

0.550000E 

03 


-0.U8756F- 

■06 

0. 500000F 

00 

0. 100000= 

01 

0.150000F 

01 

0.200000= 

01 

0.2500006 

01 

0.300000= 

01 

0.350000= 

01 

0.400000,6 

01 

0.4500006 

01 

0. 5000006 

01 

0.5500006 

03 


-0.U8756E-06 
0.500000E 00 
O.IOOOOOE 01 
0.1500006 01 
0.200000F 01 
0.250000E 01 
0.300000F 01 
0.3500006 01 
0.4000006 01 
0.4500006 01 
0.500000E 01 
0.550000F 03 


-0.118756F-06 
0.5000006 00 
0. lOOOOOE 01 
0.1500006 01 
0 . 200000 = 01 
0.2500006 01 
0.300000= 01 
0.350000= 01 
0.400000= 01 
0.4500006 01 
0.5000006 01 
0.550000E 03 


PRINTOUT OF THF MAKE CORNER COORDINATES BEFORE 
ITERATION IN THF 7.-6»I P6C.T ION 


0. 103858E-07 
0.1038906-07 
0. 103R98E-07 
0. )03890E-C7 
0. 103893E-07 
0. 1038986-07 
0. 1038986-07 
0.1038986-07 
0. 1038986-07 
0.1038986-07 
0. 1038986-07 
0. 1038906-07 


0. 1038986-07 
0.16I3898F-07 
0.103898F-07 
0.1038996-07 
0. 103898E-07 
0. 103308F-07 
0.103898C-07 
0.1038986-07 
0.1038986-07 
0.1038986-07 
0.1038986-07 
0. 103893G-07 


0. 103998F-07 
0. 1033986-07 
0. 103898F-07 
0. 103898E-07 
0. 1038996-07 
0. 103898F-07 
0. 103998=-07 
0.1038986-07 
0. 1038986-07 
0. 1038986-07 
0. 1933986-07 
0. 1038986-07 


0.103898E-07 
0.1038986-07 
0.103898E-07 
0.1038986-07 
0.1038986-07 
0.1038986-07 
0. 1038986-07 
0. 103898E-07 
0.103898E-07 
0. 1038986-07 
0.103898=-07 
0. 1038986-07 


0. 1033986-07 
0.193898F-07 
0.1038986-07 
0.103998F-07 
0. I03898E-07 
0. 1038986-07 
0. 103898P-07 
0. 1038986-07 
0.103898E-07 
0. 103898E-07 
0. 1038906-07 
0. 103398E-07 


0.1038986-07 
0.1030986-07 
0. 1038986-07 
0.1938906-97 
0. 1038986-07 
0. 10389GE-07 
0. 1D389R6-07 
0. 1O3098G-07 
0.1038986-07 
0.103898F-07 
0.1038986-07 
0.103898E-07 


0.103898F-07 
0.1038986-07 
0.1038986-07 
0, 103898=-07 
0.1038986-07 
0. 10389.86-07 
0. 103898F-07 
0.103898=-07 
0. 103898C-07 
0.10339.8F-0r 
0.1C3898=-07 
0.103898=-07 


0,1038986-07 
0.1038986-07 
0.1038986-07 
0.1038986-07 
0, I03898E-07 
0.1033986-07 
0. 103898 = -07 
0.1038986-07 
0. 1038986-07 
0.1030986-07 
0. 1038986-07 
0.103898=-07 


PRINTOUT OF THE WAK= CORMFP COOROINATFS 8=608= 
ITERATION IN TH= Y-OIR=CTION 


I 


I 


0.0 

0. 1061P2E 

01 

0. 1959186 

01 

0.26938flr 

01 

0.3265316 

01 

0.36734 76 

01 

0.3918376 

01 

0. 4000006 

01 

0.0 

0, 106U2F 

01 

0. 195918 = 

01 

0.269388= 

01 

0.326531= 

01 

0.3673476 

01 

0.3918376 

01 

0.4O0OO06 

01 

0.0 

0.106122= 

01 

0. 1959186 

01 

0.269388= 

01 

0.326531= 

01 

0.367347= 

01 

0.39)837= 

01 

0.400000= 

01 

0.0 

0. 196122= 

01 

0. 19 5913 6. 

01 

0.2693886 

01 

0.326531= 

01 

0.3673476 

01 

0.391837= 

01 

0.400000= 

01 

0,0 

0. 106 L22F 

01 

0. 195918E 

01 

0. ?693386 

01 

0.3265316 

01 

0. 36734 7F 

01 

0.301837= 

01 

0.400000= 

01 

0.0 

0, 1061226 

01 

0. 1959186 

01 

0.2693386 

01 

0. 3265316 

01 

0.3673476 

01 

0.3918376 

01 

0.4000006 

01 

0.0 

0.1061226 

01 

0. 1959186 

01 

0.269388= 

01 

0-3265316 

01 

0, 36734 7F 

01 

0.391837F 

01 

0.4000006 

01 

0-0 

0.l0612?r 

01 

0. 19 59 18F 

01 

0,269388= 

01 

0. 32.65316 

01 

0. 367347= 

01 

0.391837= 

01 

0.400000= 

01 

0. 0 

0. 106122F 

01 

0. 1959186 

01 

0.2693886 

01 

0-3265316 

01 

0-36/34 76 

01 

0- 391 837= 

01 

0.400000= 

01 

0.0 

0.106)226 

01 

0. 19 59 136 

01 

0,2693886 

01 

0.3265316 

01 

0.3673476 

01 

0.391837= 

01 

0.400000= 

01 

0.0 

0. 106 122= 

01 

0. 19591-3 = 

01 

0.2693886 

01 

0.3265316 

0) 

0.367347= 

01 

0.391837= 

01 

0.4000006 

01 

0.0 

0.1061226 

91 

0. 195918 = 

01 

0.2691886 

01 

0. 3.''65316 

01 

0.36 /34 76 

01 

0.39)837= 

OI 

0.400000= 

01 


SOI 



THf: niSTRJBUTinN OP THP OOUBLET STRFNGTH DH 
RtGHTHAND SIDE 


-0,3968'?P-02 

-0.388^9F-02 

-0.37125F-02 

-0.34288P-02 

-0.30008E-02 

-0.23927E-02 

-0.15925F-02 

-0.68218E-02 

-0.68766F-02 

-0.63782F-02 

-0.58863Et02 

-0.51416E-02 

-0.40764E-02 

-0.26464E-02 

-0.97730F-02 

-0.95619E-02 

-0.91271E-02 

-0.8't065E-02 

-0.73r)69E~02 

-0.57133E-02 

-0.35749E-02 

-O.12S60F-01 

-0. 12283E-01 

-0. U710E-01 

-0. 10753E-01 

-0.92783E-02 

-0.71273E^Q2 

-0.43533F-02 

-0. 15Q61E-01 

-0. IA720F-01 

-0,lA011E-0l 

-0.12818E-01 

-0.10967E401 

-0,82876£-02 

-0.49955E-02 

-g. 17136E-01 

-g.l6737F-0l 

' -O.lSgO'fP-Ol 

-0. 14499E-01 

tO. 12303E-01 

-0.91863E-02 

-0.54982E-02 

-0. 185539-01 

-0. 18116F-01 

-0.17189E-01 

-0.15618E-01 

-0. 13183E-01 

-0.9778 lE-02 

-0. 583476-02 
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THIS IS THE X-WAKI: VItLOCITY 


O 

kh na 

gB 


r-s 
»Tt taj 

.>K. •* 


8 


_ t-9 
^ W 
t?d 


-0- A80A06F-02 
-0- l6082Te-02 
-0-6530T0e-03 
-0. 53201AF-03 
-0,3613736-03 
-0, 258570r-03 
-0. 1916 736-03 
-0. U5832F-03 
-0.1131916-03 
-0.80A683E-04 
-0.717AA2E-04 


-0.6774096-0? 
-0.160449F-0? 
-0.847636C-O3 
-0. 52441 5E-03 
-0.3530076-03 
-0.250580P-03 
-0. 1646656-03 
-0.1399996-03 
-0, 1034616-03 
-0.C57U3F-04 
-0.6na077F-04 


-0.466116P-02 
-0.156455F-02 
-0.821966E-03 
-0.4956056-03 
-0.325844F-03 
-0, 227359 F-03 
-0. 1658226-03 
-O.125161F-03 
-0.9693106-04 
-0, 767820 F-04 
-0.6 19067P-04 


-0.4542476-02 
-0.1520606-02 
-0. 7484436-03 
-0.4322216-03 
-0.277185F-03 
-0,1914356-03 
-0.1394456-03 
-0.1057486-03 
-0.B24915=-04 
-0.6594646-04 
-0.5368916-04 


-0.4412796-02 
-0.132S5BF-02 
-0.59709ai=-03 
-0.3360616-03 
-0.2164656-03 
-0.15 176 lc-03 
-0.1126076-03 
-0.868635F-04 
-0. 6889006-04 
-0,5592906-04 
-0.4617586-04 


-0.449 1496-02 
-0.903000F-03 
-0.4032986-03 
-0.240564£_O3 
-0.1632726-03 
-0. 119269P-03 
-0.91325 16-04 
-0.7235066-04 
-0.5857016-04 
-0.4836916-04 
-0.405 120F-04 


-0.6410416-02 
-0.440632^-03 
-0. 263797f-03 
-0.1789856-03 
-0.130310F-03 
-0. 99'.833P-04 
-0.785680P-Q4 
-0,6357386-04 
-0. 5241716-04 
-0.4384776-04 
-0. 371014P-04 


0.3804516-02 
-0.2761406-03 
-0.218201F-03 
- 0 . 1589066-03 
-0.1 195606-0'3 
-0,9300896-04 
-0.743386=-04 
-0.6066986-04 
-0.5039596-04 
-0.4235706-04.' 
-0.359719F-04 


PRINTOUT OF TK6 WAKF y-VFLOCITY 


0.0 

0.0 

0.0 

0.0 

;0.0 

0.0 

0.0 

0.0 

0.0 

-0.'7 160226-04 

-0, 175023F-03 

-0.3554306-03 

-0.6431186-03 

-0.9323976-03 

-0.9948396-03 

-0.9611386-03 

0.0 

-0. 62472 lF-04 

-n. 146541P-03 

-0.264779P-03 

-0.38n6-i7P-03 

-0-45’l965F-O3 

-0.450073F-03 

-0.44O730C-03 

0,0 

-0.5490746-04 

-0.1192146-03 

-0. 190576F-03 

-0.2451476-03 

-0.264493P-03 

-0. 261428F-03 

-0.2573596-03 

0.0 

-0.4635456-04 

-0.939147F-04 

-0.1373616-03 

-0, 164443F-03 

-0.172431F-03 

-0. 170733F-03 

-0. 169176'=-03 

0,0 

-0.379370F-04 

-0.72B592F-04 

-0.100491F-03 

-0,1156616-03 

-0. 119003F-03 

-0.U8961F-03 

-0. 1 181966-03 

0. 0 

-0.305000F-04 

-0.563505F-04 

-0.7430186-04 

-0. 8463076-04 

-0. 8700276-04 

-0.8663506-04 

-0.865576F-04 

0.0 

-0.2428r7F-04 

-0.437568P-04 

-0.5675396-04 

-0.43C960E-04 

-0.652968F-04 '' 

-0.65f>978E-04 

-0.654446F-04 

0.0 

-0.1911436-04 

-0.3378236-04 

-0.430350P-04 

-0.4719376-04 

-0.479207F-04 

-0. 4755016-04 

-0.475683=-04 

0.0 

• -0.151972F-04 

-0.2660606-04 

-0.3371196-04 

-0.3704686-04 

-0. 378746F-04 

-0,3793005-04 

-0.380060F-04 

0.0 

-0.121347F-04 

-0.2114756-04 

-0.2677666-04 

-0.2958416-04 

-0.305684F-04 

•^0,3082136-04 

-0.3090566-04 

0,0 

0,0 

0,0 

0.0 

0.0 

0, 0 

0. 0 

0.0 


PRINTOUT OF the wake 7.-V6L0CITY 


-0,8748856-01 
-0.432382F-01 
-0.3530176-01 
^0.3126206-01 
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• 0.0 
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0.1186486 00 
0.119454= 00 
0.1200426 00 
0. 120485F 00 
0.1208296 00 
0.121102F 00 
0.1213^26 00 
0.0 
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